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ABSTRACT 



The purpose of this research is to design, simulate and 
implement a robust autopilot system for the vertical mode of 
operation of the Archytas prototype. Archytas is an Unmanned 
Air Vehicle that is designed to take off and land vertically, 
and to transition to horizontal forward flight. A feedback 
control scheme is designed for both the single-input, single- 
output and the multi-input, multi-output subsystems using 
optimal control techniques. In this research, the linear 
quadratic regulator performance measure is modified to allow 
for its application to the tracking problem solution. 
Additionally, the control systems are designed using reduced 
order models. Computer simulations show that the reduced 
order controller designs provide results comparable to the 
full order controller designs. Successful hardware tests with 
roll rate control system validated the reduced order model 
design philosophy used in this research. 
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I. 



INTRODUCTION 



A. THE ARCHYTAS CONCEPT 

The current inventory of Unmanned Air Vehicles (UAVs) is 
not able to meet the expanding need for real-time intelligence 
at the Marine expeditionary unit or Army battalion level. The 
Naval Postgraduate School UAV Flight Research Lab is directing 
efforts at developing a ducted-fan vertical takeoff and 
landing (VTOL) vehicle to meet these increasing needs. The 
NPS air vehicle, named Archytas, is serving as a technology 
demonstrator to evaluate the concept of a winged ducted-fan 
VTOL aircraft. The research is being directed at applying the 
technology and equipment developed in the U.S. Marine Corps' 
Airborne Remotely Operated Device (AROD) program and the U.S. 
Army's AQUILA program. 

Archytas, pictured in Figure 1.1, is designed to take off 
and land vertically. After climbing to altitude, Archytas 
will transition to horizontal flight by pitching about its 
center of gravity to a wings level attitude. The positioning 
of the duct and wings (including the canard) allow for the 
vertical takeoff and landing capability. The ability to 
transition to horizontal flight will extend the vehicle's 
range and provide the capability for a high speed dash. 
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SIP E V IEW 




Figure l.l sketch of Archytas 
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B. THE CONTROL PROBLEM 



Archytas is powered by a vertically mouirced 28-horsepower 
engine turning a three-bladed propeller. The use of a single 
propeller in a duct (ducted fan) simplifies the design, but 
intensifies the stability and control problems. The dynamic 
behavior about a given axis is coupled with other vehicle 
dynamics. In particular, there are three types of coupling: 

1. The single propeller design introduces a gyroscopic 
coupling between the pitch and yaw axes. 

2. Reactive torques are applied to the roll axis as the 
engine speed is varied. 

3. A loss of lift due to thrust occurs when the vehicle is 
pitched during the translation to horizontal flight. 

It was the goal of this thesis to design a control system 
which would allow stable flight of Archytas during takeoff and 
landing, and during the vertical mode of operation. Because 
of the coupled nature of the Archytas control problem, linear 
quadratic regulator control theory was used. 

C. THESIS ORGANIZATION 

In Chapter II, the general nonlinear equations of motion 
are developed. These equations of motion are then applied to 
Archytas with special attention given to the effects of the 
spinning propeller. Additionally, the equations of motion are 
linearized about the hover operating condition using the 
small-disturbance theory. 
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Chapter III discusses the state space representation of a 
general system and develops a procedure for selecting state 
variables for the optimal control tracking problem. Key to 
the design of the Archytas control system is the proper 
selection of the state variables. 

In Chapter IV, three control law designs are formulated 
based on the linearized hover equations. One design is for 
the single input, single output (SISO) roll rate controller. 
The second design is for the SISO altitude rate controller. 
These two SISO systems are similar in their development. 
Next, the multiple input, multiple output (MIMO) pitch and yaw 
angle controller is designed. Central to each control law 
design is the use of a reduced order model to simplify the 
design process and physical implementation. Finally, these 
control laws are applied to the nonlinear model for 
validation. 

Chapter V discusses the field test results of the roll 
rate controller. The roll rate controller was evaluated with 
Archytas mounted on a test stand to allow a roll about the 
longitudinal axis. In addition, conclusions based on the 
computer simulations, the field tests, and recommendations for 
future research are presented. 
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II. MODELING THE ARCHYTAS PROTOTYPE 
The purpose of this chapter is to develop a suitable 
dynamic model of the Archytas prototype. Because Archytas 
is a ducted fan device, special attention must be given to 
the significant gyroscopic contribution of its propeller. 

A. DERIVATION OF RIGID BODY EQUATIONS OF MOTION 

The rigid body equations of motion in this section are 
developed for the Archytas prototype in the following way. 
First, Archytas is regarded as a single rigid body, and the 
equations of motion are derived with respect to a set of 
body fixed axes. These equations are the general equations 
governing aerodynamic flight for all aircraft. Next, the 
changes introduced by the spinning rotor are evaluated and 
included in the equations. These are the gyroscopic effects 
due to the propeller. Finally, the development of a 
complete model is undertaken for the specific case of 
Archytas using the actual measurements and experimental data 
from the AROD prototype as first approximations. 

1. Force and Moment Equations 

The general equations of motion are developed for a 
typical aircraft in References 1 and 2 . A combination of 
the two approaches is taken here to arrive at the set of 
equations describing Archytas. The equations of motion are 
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obtained from Newton's second law, which states: The 
summation of all external forces acting on a body is equal 
to the time rate of change of the momentum of the body; and 
the summation of the external moments acting on the body is 
equal to the time rate of change of the moment of the 
momentum (angular momentum) . The time rates of change of 
linear and angular momentum are referred to an absolute or 
inertial reference frame. This absolute or inertial 
reference frame is an axis system fixed to the Earth. 

Figure 2.1 depicts both the body fixed axes and the inertial 
reference frame. [Ref. 1: p.84] 

Newton's second law can be expressed in the 
following vector equations: 

; ( 2 . 1 ) 

( 2 . 2 ) 

where F is the externally applied force, M the externally 
applied moment about the center of mass, y the velocity 
vector, and H the angular momentum vector about the center 
of mass. 

The vector equations, in scalar form, consist of 
three force equations and three moment equations. The force 
equations can be expressed as 
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( 2 . 3 ) 



F,= ^tau) ; : F,-^Amw) ; 

where F^, Fy, F^ and u, v, w are the components of the force 
and velocity along the x, y and z axes, respectively. The 
force components are composed of contributions due to the 
aerodynamic, propulsive, and gravitational forces acting on 
the aircraft, the Archytas prototype for the purpose of this 
thesis. The moment equations can be expressed in a similar 
manner: 



^= 4 "' '■ 






dt “ 



( 2 . 4 ) 



where L, M, N and Hy, are the components of the moment 
and moment of momentum along the x, y and z axes 
respectively . 

Now, considering Figure 2.2, let 6m be an element of 
mass of the Archytas prototype, v be the velocity of the 
elemental mass relative to the inertial axes, and let 6 f be 
the resultant force that acts upon it. Newton's second law 
then gives the equation of motion of 6m : 



5£=6in-^ . (2.5) 

dt 

The total force acting on the vehicle is a summation of all 
the forces that act upon all the elements. The internal 
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Figure 2.2 An Element of Mass on the Archytas 
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forces, those exerted by one element upon another, all occur 
in equal and opposite pairs, by Newton's third law, and 
hence contribute nothing to the summation. Thus U6F=F is 
the resultant external force acting upon the vehicle. The 
velocity of the differential mass 6m is given by 



V=V 

^ dt 



( 2 . 6 ) 



were is the velocity of the center of mass of the 
aircraft and dr/dt is the velocity of the element relative 
to the center of mass. Substituting this expression for the 
velocity into Newton's second law, equation (2.1), yields 



. < 2 - 7 ) 

Assuming that the mass of the aircraft is constant, equation 
(2.7) can be rewritten as 



F=J-y (V +-3^) 6m 



( 2 . 8 ) 



or 



dv 

F=m—^ 

dt 



+ — - rbm . 



( 2 . 9 ) 



Because r is measured from the center of mass, the summation 
Xr.6m is equal to zero and the force equation (2.9) becomes 
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( 2 . 10 ) 



Which relates the external force on the aircraft to the 
motion of the vehicle's center of mass. 

The relationship between the external moment and the 
rotation of the aircraft is obtained from a consideration of 
the moment of momentum. For the differential element of 
mass, 6m, the moment of momentum is by definition 6H=rxvbm, 
The moment equation can be written as 



6M=-^6H=-^(rxv)6m . ( 2 . 11 ) 

dc dt 

The velocity of the mass element can be expressed in terms 
of the velocity of the center of mass and the velocity of 
the mass element relative to the center of mass: 



v=v + 

— C 



dt 






( 2 . 12 ) 



where w is the angular velocity of the vehicle and r is the 
position of the mass element measured from the center of 
mass. The total moment of momentum can be written as 



(Mxr) ] 6t? . (2.13) 



The velocity is a constant with respect to the summation 
and can be taken outside of the summation sign 
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H=(^r8/7?) X (rxo)xr) 6/n . 



(2.14) 



The first term in equation (2.14) is zero because the term 
5m=0 as explained previously. The position vectors and 
angular velocity can be expressed as 

r=xi+yj+zk ; (2.15) 

co=pi+gjL+ric ; (2.16) 

where p, q and r are the scalar components of w, and i, i, k 
are unit vectors in the directions of x, y and z. 
Substituting w and r into equation (2.14) and expanding, H 
can be written as 

H= (pi+<zi+rk)'^ {x^+y^+z^) bm-^ (xi+y^+zk) {px+gy+rz) bm . 



The scalar components of H are 



(2.17) 



H^= pY,{y^*z^) bm-qYxy bm-zY^xz bm ; 

Hy=-pY^y bm+qYi^^*^^) S/n ; (2.18) 

-pY^z bm-qYyz 6/7}+r£ (x-2+y^) bm . 

The summations in the above equations are the mass moment 
and products of inertia of the aircraft and are defined as 
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(2.19) 



^x=/ / / 5^ , f fw 5^ ; 

Iy=JJj{x^->-z^)6m, I^^=Jjjxz dm ; 

V V 

(^^+y^)6^ / Iy^=fJJyz dm . 

V V 

The terms ly and are the mass moments of inertia of 
the body about the x, y and z axes, respectively. The terms 
with the mixed indices are called the products of inertia. 
Both the moments and products of inertia depend on the shape 
of the body and the manner in which its mass is distributed. 
The larger the moments of inertia the greater the resistance 
the body will have to rotation. Applying the notation of 
equation (2.19) to equation (2.18), yields the scalar 
equations for the moment of momentum: 

^x~ P^x P^xy ^^xz ' 

Hy=-Pl^^qly-Zly, ,' ( 2 . 20 ) 

^z~ P^xz ^^yz^^^z • 

If the reference frame is not fixed to the aircraft, then, 
as the aircraft rotates, the moments and products of inertia 
will vary with time. To avoid this difficulty an axis 
system will be fixed to the aircraft (body axis system) . 

Now the derivatives of the vectors y and H referred to the 
rotating body frame of reference must be determined. 
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It can be shown that the derivatives of an arbitrary 
vector A referred to a rotating body frame having an angular 
velocity o) can be represented by the following vector 
identity 



dA 

dt 



dA 

dt 



+C0.XA ; 



( 2 . 21 ) 



where the subscripts I and B refer to the inertial and body 
fixed frames of reference, respectively. Applying this 
identity to equations (2.1) and (2.2) yields 



dy.. 

E=m-^ +ir?(a)xv^) ; 
UC g 



( 2 . 22 ) 



M= 



dH 
dt , 






( 2 . 23 ) 



These are the general equations governing aerodynamic flight 
and have the scalar components: 



F^=m(u+gw-rv) , F =m{v+ru-pw), F^=m{w+pv-qu) ; ( 2 . 24 ) 



L=H^+qH^-rHy, M=Hy+iH^-pH^, N=H^+pHy-qH^ . ( 2 . 25 ) 

The components of the force and moment acting on the 
aircraft are composed of aerodynamic, gravitational and 
propulsive contributions. 

At this point, it is recognized that most aircraft 
have a plane of symmetry. If the xz plane is selected to 
coincide with this plane of symmetry, then from equation 
(2.19), l^=ly^=0 must be satisfied. However, for the case of 
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the Archytas , the three blades of the propeller provide two 
planes of symmetry, xz and xy. Thus, the products of 
inertia, I^, and I„, equal zero. The moment equations 
for Archytas can now be written as 

M=Iyq+rp(I^-I^) ; (2.26) 

N=I^t+pq(Iy-I^) . 

2. Effect of the Spinning Rotor 

Archytas, like AROD, is a gyroscope. The single 
propeller rotates about the longitudinal vehicle axis to 
produce a downwash or jet of air through the duct which 
makes up the Archytas body. This spinning rotor exerts a 
gyroscopic moment on the body of the vehicle. Reference 3 
states that in developing the equations of motion for 
aircraft with propellers which exert gyroscopic moments on 
the body "more often than not, such gyroscopic moments turn 
out to be negligible." However, as demonstrated by Bassett 
[Ref. 4 : p. 19], in the AROD case the angular momentum of 
the propeller, Ip(ARoo)/ equals 11.3 ft^-lb^/sec . Compared 
with AROD's nominal total mass of 2.64 Ib^, (85 Ib^) , it is 
clear that the angular momentum imparted by the propeller is 
significant and that gyroscopic effects will play a large 
part in modeling the dynamic behavior of AROD. For 
Archytas, Ip is identically equal to Ip(AROD>' H-3 ft^-lb„/sec. 
Compared with Archytas' nominal total mass of 3.11 Ib^, (100 
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Ib^) , it is clear that similar to AROD the gyroscopic 
effects will be significant and must be included in modeling 
the dynamic behavior of Archytas. This gyroscopic moment 
can be accounted for as follows. 

Angular momentum, Hp, due to the propeller is 
defined as 

'■ ( 2 . 27 ) 

where Ip is the propeller moment of inertia, and cjp is the 
propeller angular velocity. Since the propeller lies in the 
yz plane and spins symmetrically about the x axis, Hp is 
directed only along x and Hpy=Hp^=0. Equation (2.27) becomes 

• < 2 . 28 ) 

Etkin [Ref. 2; p.93] states that the resultant angular 
momentum of an aircraft with spinning propellers is obtained 
simply by adding Hp to the H previously defined by equation 
(2.20). Adding equations (2.28) and (2.20) and keeping in 
mind that the products of inertia equal zero, yields 

Hy=qly ; ( 2 . 29 ) 

Applying equation (2.29) to equation (2.25), the moment 
equations can be written for the specific case of Archytas 
as 



16 



L=I^p+qriI^-Iy] +Ip(i)p ; 

Af=J^g+rp( Jjf-Jj) +rJpCOp ; (2.30) 

N=l^t+pg{ly-IJ -gJpOJp ; 

and the force equations are those of equation (2.24). 

3. Orientation and Position 

Because the frame of reference developed for the 
equations of motion is fixed to the aircraft, and moves with 
it, the position of the aircraft cannot be described 
relative to it. The orientation and position of the 
aircraft can be defined in terms of a fixed frame of 
reference as shown in Figure 2.3. [Ref. 1: p. 89] 

The orientation of the aircraft can be described by 
a series of three consecutive rotations, whose order is 
important. The angular rotations are called the Euler 
angles. The orientation of the body frame with respect to 
the fixed frame can be determined in the following manner. 
The aircraft is imagined first to be oriented so that its 
axes are parallel to the fixed frame and the following 
rotations are then applied. [Ref. 2: pp. 89-91] 

1. A rotation about oz,, carrying the axes to CX 2 Y 2 Z 2 
(bringing Cx to its final azimuth) . 

2. A rotation 0 about oy 2 , carrying the axes to CX 3 y 3 Zj 
(bringing Cx to its final elevation) . 

3 . A rotation # about 0X3, carrying the axes to their 
final position Cxyz (giving the final angle of bank to 
the wings) . 
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Figure 2.3 Relationship between Body and 
Fixed Axes System 
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Now that the Euler angles are defined, the flight velocity 
components relative to the earth-fixed reference frame can 
be determined. To do this, let the velocity components 
along the x, , y,, z, frame be dx/dt, dy/dt, dz/dt and, 
similarly, let the subscripts 2 and 3 denote the components 
along x,, Y 2 , and Xj, yj, Zj, respectively. From Figure 
2.3, it can be shown that 



where 



and 



and 



dt ^ 



dt 



dz ^ 
dt 






(2.31) 



Ui=U2COSilr-V2sini(r ; 

Vj^ = U2sin\|;+V2COS\|r ; (2.32) 



U2 = u3cos0 + W3sin0 ; 

^2 = ^3 ; (2.33) 

W2 = -U3sin0+W3cos0 ; 

U3 = U ; 

V3 = vcos^»-sin<I> ; (2.34) 

W3 = vsin^>+wcos^> ; 



from this, the absolute velocity in terms of the Euler 
angles and velocity components in the body frame can be 
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determined. Note the shorthand notation S^ssin'I^, C,^=cos^, 
S 0 =sin 0 , etc, used in the following equations: 







dt 




dy 




dt 




dz 











S^SqC^-C^S^ 
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( 2 . 35 ) 



Integration of these equations yields the aircraft's 
position relative to the fixed frame of reference. The 
relationship between the angular velocities in the body 
frame (p, q and r) and the Euler rates ( ^ and ) can 

also be determined from Figure 2.3. 



( 2 . 36 ) 



Equations (2.36) can be solved for the Euler rates in terms 
of the body angular velocities and is given by equation 

(2 . 37) 
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( 2 . 37 ) 



0 S^secQ C^secd 

By integrating the above equations, the Euler angles (0, ij/ 
and 0 ) can be determined. 
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4. Gravitational and Thrust Forces 



The gravitational force acting on the aircraft acts 
through the center of gravity of the aircraft. Because the 
body axis system is fixed to the center of gravity, the 
gravitational force will not produce any moments. However, 
the gravitational force will contribute to the external 
force acting on the aircraft and will have components along 
the respective body axes. From Figure 2.4 the gravitational 
force component in the direction of each axis is found to be 

X^=-mg cosGcosT ; 

Yg= /ng cosGsinT ; (2.38) 

Zg=-mg sin0 . 

With the aerodynamic forces (including the propulsive 
forces) denoted by (X, Y, Z) , the resultant external forces 
are 



F^=X-mg cosGcosT; 

Fy=Y+mg cos6sin}P ; (2.39) 

F^=Z-ing sin6 . 

5. Summary of Equations of Motion 

In the previous sections, the equations that 
completely describe the dynamic behavior of Archytas have 
been developed. Equations (2.24) and (2.30) define the 
externally applied forces and moments which are represented 
by F^, Fy, Fj and L, M, N. Through the Euler angles, the 
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behavior of Archytas can be observed relative to the Earth. 
Specifically, the translational velocities for the fixed 
frame of reference, dx/dt, dy/dt and dz/dt, can be 
determined from the body-fixed velocities, u,v,w and the 
Euler angles, 0 , ^ and $, using the transformation of 
equation (2.35). Additionally, equation (2.37) describes 
the relationship between the Euler angles and the body 
angular velocities, p, q and r. Table 2.1 gives a summary 
of the rigid body equations of motion. 

B. ARCHYTAS NONLINEAR SYSTEM EQUATIONS 
1. Applied Forces and Moments 

The Archytas model is developed using the 
measurements and data from the AROD prototype as first 
approximations. These measurements and data were obtained 
by the AROD project engineers based on wind tunnel tests. 
This data consists of tabular results that describe the 
aerodynamic lift and drag coefficients and physical 
measurements of constants such as weight, moments of inertia 
and servo gains. This tabulated data forms the basis from 
which the applied forces and moments may be determined. The 
forces and moments, which are computed from the data are of 
two types: aerodynamic and thrust. This data is listed in 
Appendix A. 
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TABLE 2.1 SUMMARY OF GENERAL EQUATIONS OF MOTION 



X-mg = m((i+qw-rv) 
Y+mg CgSjf = m(v+zu-pw) 
Z-mg Sq - m(w+pv-qu) 

L=IJ>+qr{I^-Iy) +IpCip 
M=Iyq+ip(I^-I^) *rIpWp 
N=IJ:-^pq ( ly- 1^) - gJpUp 

p=4*-tSe 

g=ec*+TCe5* 

r=^C^C^-dS^ 

0=gC^-rS^ 

^=p+qS^TQ+rC^T^ 

T= (gS^+rCj) sec0 



Force equations 



Moment equations 



Angular velocities 



Euler rates 



Velocity of aircraft in the fixed frame in terms of Euler 
angles and body velocity components 



dx 




dt 




dy 




dt 




dz 




dt. 
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Sq 




SqS^ + 











V 



w 



a. Total Angle of Attack and Body Roll Angle 

The aerodynamic data describes the forces and 
moments relative to the vehicle's total velocity vector, 

Vtot* These forces and moments (in the body-fixed coordinate 
system) depend on the total angle of attack {a') and the 
body roll angle (A) . The total angle of attack and body 
roll angle can be defined in terms of the velocity 
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components as shown in Figure 2.5. The equations for a' and 
A are given as: 



a 



=sin 



-1 



vwterm 

u 



(2.40) 



and 



A=tan"^— . (2.41) 

w 

b. Aerodynamic Forces and Moments 

The forces computed from the tabular data are 
lift (F,) and drag (Fj) . These forces are depicted in Figure 
2.6(a). The transformation from lift and drag to F^^, F^y and 
F„ is given as 



and 
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-cosa ' 



-cosa 
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(2.42) 







sin A 






cosA_ 



ayz 



(2.43) 



where F„, F^y and F,^ are the forces in the x, y and z 
directions due to the aerodynamic data. 

Similarly, the aerodynamic moments applied to the 
body axes as a result of the vehicle's movement through the 
air can be derived. These moments (shown in Figure 2.6(b)) 
are referred to as the aerodynamic angular moments of roll 
(Lj) / pitch (M,) and yaw (N,) and are given as 
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y^-z^ plane 




Figure 2.5 Angle of Attack (a' 
Body Roll Angle (A) 



► X, 



) and 
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(b) Moments 

Figure 2.6 Aerodynamic Forces and Moments 
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(2.44) 



' La ' 




sina " 


-cosa ' 




M.; 


^pyyz. 




-cosa " 


-sina 




— 1 



and 







sinA 






cosA. 



where and are the yaw and roll moments relative to 
VjoT* Mpyyj, is the moment in the y-z plane. The relationships 
that define F,, and My, M, were developed by the AROD 

project engineers and are listed in Appendix A. [Ref. 15] 

c. Control Forces and Moments Due to Command Inputs 

The forces and moments previously discussed are a 
result of the vehicles motion through the air. A second 
category includes those forces and moments which are a 
result of the commanded inputs. These inputs control the 
rotor speed and the displacement of the control surfaces. 
[Ref. 4: p. 33] 

(1) Forces Due to Induced Thrust of the Ducted 
Fan. The commanded inputs include the ability to change the 
rotor speed; thus, changing the thrust provided by the 
ducted fan. Due to the orientation of the body-fixed axes, 
the force due to the thrust (F-n,,) is directed completely 
along the x-axis. The relationship that defines F-n,, was 
developed by the AROD project engineers and is given as 
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( 2 . 46 ) 



F^t^^ = THOVER + XRPM * 6 .^^ 

where THOVER is a constant equal to the nominal thrust at 
hover. THOVER is set equal to Archytas prototype weight of 
76.5 ft/ lbs. XRPM is the slope of the thrust versus engine 
rpm curve. For the hover rpm of 712.0943 rad/ sec (6800 
rpm) , XRPM is equal to 0.2387 Ibf/rad/sec. h.p^ is the 

change in engine rpm from the nominal hover rpm in rad/ sec. 

(2) Moment Due to Ducted Fan Effects. A 
gyroscope imparts no torque on its axis if it spins with a 
constant angular rotation. In hover (constant rotor speed) , 
Archytas behaves similarly to a gyroscope. If the rotor 
accelerates (positively or negatively) a torque is applied 
to the axis. This torque is accounted for in Equation 
(2.30). However, Archytas is a ducted fan and the drag 
between the rotor tip and the inside body wall creates a 
moment about the x-axis (roll) . The project engineers for 
AROD determined an approximation for this moment based on 
experimentation. Because the Archytas duct is identical to 
the AROD duct, the moment determined for AROD applies to 
Archytas and is given as 

= ^..cc^rpn, < 2 . 47 ) 

where 5^^,^ is defined above. is a constant which is 

dependent on the duct geometry and is equal to 0.0729. [Ref. 
4; pp. 33-34] 
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(3) Moments Due to Control Surface Displacements . 
The commanded inputs also include the displacement of the 
four control surfaces (vanes) within the downwash from the 
duct. The displacement of these vanes within the downwash 
of the duct imparts moments about the body axes. Figure 2.7 
shows how the vanes are symmetrically arranged below the 
duct. The vanes are displaced by a servo mechanism 
connected directly to the top of each vane. Vanes (1) and 
(3) are operated together as "elevators" and impart a moment 
about the y-axis (pitch). Vanes (2) and (4) together are 
the "rudder" and contribute a moment about the z-axis (yaw) . 
Vanes (1) and (3) displaced in opposite directions and (2) 
and (4) displaced oppositely work as "ailerons" to impart a 
moment about the x-axis (roll) . The actual torque applied 
by each combination of vanes was determined experimentally 
and described by "constants of effectiveness" which were 
calculated by the AROD project engineers. These constants 
of effectiveness can be applied directly to Archytas . [Ref. 

4 : p. 34 ] 

The constants of effectiveness are given the symbols 
Ljjff, M„ff and N„ff, for their contribution of moments about the 
roll, pitch and yaw axes due to a displacement by the 
ailerons, elevator and rudder. The relationships resulting 
in moments about the three body axes are 
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Figure 2.7 Archytas control Vanes 
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( 2 . 48 ) 



r = r X 

M -M -- P. S • 

^^eerr ^ factoz^ e ' 

^t=^reff^faccor^r > 

where 6^, 6g and 6^ are the displacements of the aileron, 

elevator and rudder respectively. and are scaling 

factors. The relationships that define and were 

developed by the AROD project engineers and are listed in 
Appendix A. K^ff and N„ff equal -150,379.57, -112,716.87 

and -128,774.80, respectively, with units of lbn,-in‘- 
rad/sec^. Because of the symmetry of Archytas, cross 
coupling of the control surfaces is negligible and is 
ignored. [Ref. 4: pp. 35-36] 

2. Servo Equations for Control Surfaces and Throttle 
The model airplane servos used to drive the control 
vanes and throttle linkage are identical and can be modeled 
as second-order dynamical systems. The response of these 
servos to a step response was measured. These measurements 
were used to compute the natural frequency and damping 
ratio. The results are summarized below in the form of 
constants H, and Hj. A detailed explanation as to how these 
results were obtained is contained in Appendix B. 
a. Control Surface Servos 

Three equations will describe the operation of 
the elevators, rudders and ailerons. For each of these 
equations, at least two servos are operating at the same 
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time on different control vanes. Each servo receives a 



command input (pulse-width, modulated signal) which results 
in an angular displacement of the servo. The corresponding 
differential equations for the servos can be written as 

; ( 2 . 49 ) 

where H, and H 2 are the servo gain constants equal to 71.1 
and 2745.8. u,, u, and u^ are the servo inputs. 6^, 6^ and 

6j. are the servo position angles. 

b. Throttle Servos and Engine 

The servo motor used to open and close the 
throttle is identical to the servo motors used for the 
control surfaces. The 2-cycle, 2-cylinder gasoline engine 
with dual carburetors can be modeled as a first order lag 
system. The complete third order system is 






( 2 . 50 ) 



where Wg is the lag time constant equal to 2.0 rad/sec and 
Kg is a scaling factor equal to 837.758 rad/sec/rad. 

3 . Summary 

The result of the previous section was nine 
equations completely describing the dynamic behavior of 
Archytas. These equations, combined with the equations that 
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define the servos and engine, form the complete nonlinear 
system. This nonlinear system will be the base model 
throughout this thesis. Table 2.2 lists the equations 
rearranged so that the dynamic variables of interest may be 
determined. Additionally, the applied forces and moments 
defined in this section have been substituted into the 
equations. These applied forces and moments complete the 
development of the model for the specific case of Archytas. 

C. LINEARIZING THE TJICHYTAS MODEL 

There exist many analytical and graphical techniques for 
controller design and analysis of linear systems. 

Conversely, there are no good methods available for solving 
a wide class of nonlinear systems. Thus, in the design of 
control systems it is practical to first design the 
controller based on the linear system model generated by 
neglecting the nonlinearities of the system. These 
nonlinearities are neglected by linearizing the model about 
a steady-state reference condition. The designed controller 
is then applied to the nonlinear system model for validation 
and subsequent redesign if necessary. In this section, a 
linear model is generated from the equations of Table 2.2 
based on steady-state assumptions and physical 
approximations. [Ref. 5: p. 11] 

The nonlinearities of the Archytas system equations fall 
into two categories: (1) nonlinear combination of states 
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Table 2.2 ARCHYTAS NONLINEAR SYSTEM EQUATIONS 



[L^ (6_.pJ +L, (63) +{Iy-l 2 ) rg-J.^(bp] 

■“X 

^=4’ (5^) + ( J,- Jj pr- J.^Wpr] 

*^7 

f=A +iv,(5j + (J^-Jj,)pg+J^^0)pg] 

*^2 

ii= (rv-grv^) -gcos0cosT + A F^) +F 7 .^rx(^rpai) 3 

v={pw-ru) +gcos0sinT + — F,„(F., FJ 

777 

w=(qu-pv) -gsin0+— F 3 , (F,, F.) 

m 

0 =gcos$-rsin^> 

^=p+gsin4>tan0+rcos^»tan0 
'P= (gsin^»+rcos^») sec 0 

t~H2^t^H2Ut 

^ rpa! ^£^rpin''’^£^F^t 

(e.g. Equation (2.30)) and (2) discontinuous functions (e.g. 
table lookup of aerodynamic force and moment coefficients) . 
These nonlinearities will be neglected and the nonlinear 
equations will be replaced with linear approximations. 

1. Steady-State Assvunptions 

The motion of Archytas in the hover mode consists of 
small perturbations from a steady-state condition. The 
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steady-state hover condition is defined as a’ equal to 90°. 
All translational and angular movement is very small. The 
steady-state hover condition results in the following 
simplifications of the nonlinear system equations. 

1. The product of two small numbers is an extremely small 
number, thus terms involving the products of 
translational or angular velocities are equal to zero 
(e.g. rq, pq, (pw-ru) are set equal to zero) . 

2. The aerodynamic forces (F,^, F,y, F„) and moments (L^, M,, 
Nj) are very small for a' equal to 90°, and can be 
neglected in the hover flight condition. (Note: 
Because the vanes lie within the downwash of the 
propeller, the vane effectiveness coefficients can not 
be neglected.) 

3. The sine of a state is equal to the state and the 
cosine of a state is equal to one. This is the small 
angle approximation for angles less than 15 degrees. 

4. Kjuc, is a small number equal to 0.0729. For hover or 

very small translational velocities, is a very 

small number. Therefore, the product of K<,uct ^rid 6^p„, 
is neglected. 

5. The propeller angular velocity, (Op, is considered 
constant. Thus, the propeller angular acceleration, 
(bp, is equal to zero. 



Table 2.3 lists the Archytas system equations when the above 
simplifications are applied to the nonlinear system 
equations. Note that much of the coupling between states 
and all of the nonlinear products of states have been 
eliminated. 
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Table 2.3 ARCHYTAS LINEARIZED HOVER EQUATIONS 



P=L,(6J 

r=N,(b,) 



v= gW 
w=-gd 

6=g 

4>=p 

'f=r 

h,=-H^t,-H2b,^H2U, 

6,=-H^6,-H2b,^H2U, 

^ rpm~ ^ zpm^^ t 



2. Physical Approximations of Force and Moments 

The force, Fy^rx/ in Table 2.3 is a function of the 
engine rpm and is computed by equation (2.46). The moment 
terms, L^, M, and N,, are functions of the displacement of the 
aileron, elevator and rudder. They are computed by equation 
(2.48). The moments of inertia, I^, ly, and I„, were 
determined by the AROD project engineers and are listed in 
Appendix A. The angular velocity of the propeller, cOp, in 
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the hover mode equals 712.0943 rad/sec (6800 rev/min) . The 
v/eight of the AROD prototype is equal to 7 6.5 lbs. Using 
these constant parameters, the forces and moments for the 
hover steady state reference condition were computed and are 
listed in Table 2.4. Substituting the results of Table 2.4 
into the equations of Table 2.3 yields the linearized hover 
equations, which are summarized in Table 2.5. 

D, SUMMARY 

The nonlinear differential equations of motion of a 
rigid body were developed from Newton's second law. The 
equations were then modified for the specific case of 
Archytas. This modification included the effect of the 
spinning three-bladed rotor. The equations were linearized 
about a steady-state hover condition. Next, the applied 
forces and moments, and the servo and engine equations were 
defined. The forces and moments specific for the hover 
reference condition were computed and applied to the 
Archytas linear model. 
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Table 2.4 FORCES, MOMENTS AND CONSTANTS 



Momenta 

^t=^aeff^a 

L^=-150,379 . 5763 

~ ^eeff^fa ccor^e 

(-112,716 . 87) ( 1 ) 6 ^ 

~^zeff^ factor^ r 

^^(.= (- 128 , 774 , 80 ) ( 1 ) 6 ^ 

Moments of Inertia 
J^=7063.39, Jy=7768.22, J^=7729.58, J,^=69.552 [Ib^ in^) 

Angular velocity 

(Si =112 . 0943^^; 

^ sec 

Weight /gravity /mass 

weight=76.5 lbs.; giavi ty=3 2 . 17 4 ; mass= 

sec^ g^^vizy 

Engine Lag Time Constant and Scal e Factor 

0 )^= 2. 0 Ke= 837 .758 rad/ sec/ rad 

sec 



Forces 

F.^^^^=THOVER + XRPM 
F.^^=7 6.5+0.23 87 6_,,„ 
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Table 2.5 LINEARIZED HOVER STATE EQUATIONS 



iTATE 



-150,379.57 s _ 
^ 7063.39 



g= 



r= 



7768 . 22 
1 

7729 . 58 



-14 . 516^-6 



[-112,716 . 876^-52,440 . 82r] 

[-128,77 4 . 806, +52, 440 . 82g] =-16 .6 86, +6 



. (76.5+0.23876,„J 

“ = 2738 ^- 32 . 174 = 0 . 100296 ^^^ 

V= 32.174T 
v=-32 . 1740 



0=g 

^=p 

'f=r 

5^=-h,6,-//26,+//2U, 

5 t~^2^t~^^2^t 

^ rpm ^ Tpm^^ t 



.75r 

.78g 
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III. OPTIMAL CONTROL THEORY 

The purpose of this chapter is to develop a procedure 
for applying optimal control theory to the solution of 
tracking problems. First, several reasons for desiring an 
optimal solution are presented. Next, the state space 
representation (both continuous and discrete) of a system is 
given. Finally, the application of optimal control to the 
solution of tracking problems is illustrated. 

A. WHY OPTIMUM CONTROL FOR ARCHYTAS? 

The first reason for seeking an optimum controller is 
that feedback gains can be computed for a much broader range 
of control problems. Specifically, optimal control provides 
solutions for high order, multiple-input, multiple-output 
(MIMO) systems. Such systems are often intractable with 
classical methods. The pitch and yaw angle controller for 
Archytas is an eight order MIMO system. Thus, optimal 
control is the preferred method. 

Additionally, optimal control lends itself nicely to a 
discrete time solution of the control problem. Archytas 

employ an on board digital computer to perform inflight 
stability and control. While a continuous time controller 
can be easily discretized in many cases, design of a sampled 
data controller will simplify the procedure. [Ref . 3: p. 58] 
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Finally, optimal control is known to provide robust and 
insensitive solutions to the feedback control problem. 
Assuming that an appropriate performance measure is chosen 
to determine the optimal feedback gain matrix, K, the 
solution can be expected to have a fair degree of tolerance 
to plant model inaccuracies. Clearly, robustness is not 
only a desired property of the controller, it is an absolute 
necessity if the controller is to be applicable to both the 
linear and nonlinear models of Archytas . [Ref . 3: pp. 58-59] 

B. STATE SPACE REPRESENTATION 

The state of a system may be defined as the minimum 
amount of information required such that (given the input to 
the system) the response of the system is completely 
determined for all future time. For dynamic systems, the 
response of the system is defined by the differential 
equations that model the system, the initial conditions and 
the forcing function. The number of state variables or 
states is equal to the total order of the systems 
differential equations. In order to provide a systematic 
mathematical approach to analysis of the system, it is 
convenient to describe the system by a set of first-order 
differential equations. This set of equations is called the 
state equations and constitute the basis for the state space 
representation. [Ref. 8: pp. 206-207] 
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1. Continuous Time Systems 

The state space representation of a general n* order 
continuous-time, time-invariant system is described by the 
following matrix state equations; 

x{t) =A Kit) +M u{t) ; (3.1) 

y:{t)=Qx{t) ; (3.2) 

where the definitions in Table 3.1 apply to a system with I 
control inputs and m measurable outputs. The equations 

TABLE 3.1 



STATE SPACE DEFINITIONS FOR CONTINUOUS -TIME SYSTEMS 



Term 


Dimension 


Definition 


Ki t) 


in X 1) 


State vector 


Xit) 


[m X 1) 


Output vector 


4 


[n X n) 


Plant matrix 


s 


in X H) 


Control distribution matrix 


£ 


(m X n) 


Output Distribution matrix 



listed in Table 2.5 are linear and time-invariant. 

2 . Discrete Time Systems 

As was noted earlier, there are many benefits for 
seeking a discrete time solution for the Archytas control 
problem. Therefore, the automatic control systems designed 
will focus on the application of optimal control theory to 
discrete time systems. 
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similar to the continuous-time system, the state 
space representation of a general n‘^ order discrete-time 
system is described by the following matrix state equations: 

x{k+l) x{k) +'£ u{k) ; (3.3) 

Zik)=£x(k) ; (3.4) 

^ and C are defined as: 



; (3*5) 

r=fV^dtS; (3.6) 

~ Jo . ~ 

where T is the sampling period and k is an integer time 
index. Reference 10 provides a detailed development of the 
relationship between continuous-time and discrete-time 
systems. The definition in Table 3.2 apply to a system with 
i. control inputs and m measurable outputs. 

TABLE 3.2 



STATE SPACE DEFINITIONS FOR DISCRETE-TIME SYSTEMS 



Term 


Dimension 


Definition 


2 c(ic) 


{n X 1) 


State vector 


}:(k) 


(m X 1) 


Output vector (0 < m < n) 


£ 


(n X n) 


Plant matrix 


r 


(n X {) 


Control distribution matrix 




(m X n) 


Output Distribution matrix 
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C. LINEAR QUADRATIC REGULATOR PROCEDURE 

The theory behind the' linear quadratic regulator is well 
developed. Procedures exist which make the design of 
controllers using linear quadratic regulator theory easily 
obtainable. The purpose of this section is to illustrate 
the application of the linear quadratic regulator to the 
tracking problems associated with Archytas. Reference 11 
provides a detailed development of the linear quadratic 
regulator theory. 

1. Quadratic Cost Function 

The discrete LQR synthesis problem is that of 
determining the control that minimizes the performance 
measure: 



J=^ x'^ik) Q x(k) + u'^ik) R u{k) ; (3.7) 

k=o 



where 



J = Scalar cost of operating the system; 

x(k) = State vector at discrete times; 

u(k) = Control vector at discrete times; 

Q = Symmetric state trajectory weighting 
matrix; 

R = Symmetric control weighting matrix; 

T = Matrix transpose operator. 
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The solution to this optimization problem is the linear 
controller: 



u{k) =- ^ x{k) - -K K^k) ; (3.8) 

where M satisfies the algebraic Riccati equation (ARE) : 

g = M A + - M E + n ■ (3.9) 

[Ref. 11: pp. 345-346] Many software packages, including 
the program MATLAB used in this thesis, contain subroutines 
that calculate the value of K for a given dynamic system and 
performance measure. 

2. Performance Weighting Matrices 

The optimal control is fully specified by the system 
model and the weighting matrices Q and R. Q penalizes 
deviation of the state vector x from the origin and R 
penalizes the use of too much control effort. Trial and 
error was used in selecting values for these performance 
weighting matrices. The relationship between the weighting 
matrices Q and R and the dynamic behavior of the closed-loop 
system depends of course on the matrices A and B and is 
quite complex. The approach taken in the design of the 
controllers for Archytas was to solve for the gain matrices 
K that result from a range of weighting matrices 2 and R, 
and then simulate the corresponding closed-loop response. 

The gain matrix K that produced the response closest to 
the desired design objectives was chosen. [Ref. 11: p. 341] 
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3 . Optimal Tracking Systems 

The regulator and tracking problems appear very 
similar. The tracking problem attempts to drive the states 
of the system to a desired level; whereas, the regulator 
attempts to drive all of the states to zero. Their 
differences present conceptual difficulties and sensitivity 
problems when viewed in a practical context [Ref. 12: pp. 
643-647]. Burl [Ref. 13] presents a comprehensive 
development of the subtleties encountered when applying LQR 
synthesis to the tracking problem. This development is 
generalized below to demonstrate these subtleties and to 
indicate a procedure for selecting the proper form of the 
performance measure. The following development is applied 
to a first order system, but the results are applicable to 
systems of arbitrary order. 

Given the scalar plant: 

y ( t) =-Ay ( t) +Su( t) . (3.10) 

The purpose of the control system is to drive the output 
y(t) to the constant reference input r. This results in the 
error equation: 

e( t) =r-y ( t) ; (3 • H) 

where the desire is to drive e(t) to zero. The application 
of linear regulator theory requires that this error be a 
linear combination of the states of the system. This can be 
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readily accomplished by generating a new state equation for 
e(t) : 

e{C) =-y(t) =Ay{C) -Bu(t) ; (3.12) 

adding and subtracting Ar yields: 

e{C) =Ay(C) +Ar-Ar-Bu{t) ; 
e(t) =A[-z+y{t) ] -Bu(t) +Ai ; 

e ( t) = -Ae ( t) -Bu ( t) +Ar . (3.13) 

Since the error equation, Equation (3.19), is linear and 
time invariant, the performance measure: 

J=f{Qe^ it) it)} dt ; (3.14) 

0 

should provide an optimal solution for the system of 
Equation (3.21). However, as shown by Burl [Ref. 13] this 
optimal control problem will not yield a solution. To 
demonstrate this fact, note that the existence of this 
integral requires that both u(t) and e(t) approach zero as t 
tends to infinity. If e(t) approaches zero, then e(t) must 
also approach zero which from Equation (3.21) yields: 

0 = -A(0) -Buie) -* U(t) =4^ 0 (3.15) 

B 

This nonzero value for u(t) implies that the performance 
measure is infinite. A solution to this problem would be to 
formulate the performance measure as 
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‘J'=f lC)e^(C) + lu( t) dc 



(3.16) 



This will result in the existence of a theoretically optimal 
control provided that the model of the system (the 
parameters A and B) is exactly known. When errors exist in 
the model, a nonzero steady state error will exist. The 
resulting control system will be mathematically correct, but 
it will be unacceptable for many applications due to its 
sensitivity to changes in the plant parameters. The control 
system will be super-tuned (it will not be robust) . [Ref . 12: 
p. 645] 

This undesirable result can be overcome by letting the 
controller find the appropriate steady state value of the 
control. This is achieved by application of the performance 
measure: 



The control that minimizes Equation (3.25) given the system 
of Equation (3.21) is found by first differentiating 
Equation (3.21), yielding: 



Noting that the error is the integral of e{t) results in 
the state space system: 




(3.17) 



0 



e( t) =-Ae( t) -Bu( t) 



(3.18) 
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This system and the performance measure of Equation (3.25) 
form a linear regulator problem whose solution is state 
feedback: 



u{t) =- [k,_ 



k,] 



e{ t) 
e( t) 



(3.20) 



The gains k, and kj are computed by application of linear 
regulator theory. The actual control that is applied to the 
system is found by integrating Equation (3.28); 



u(C) =-k^je{T:) dx-k2e(C) (3.21) 

0 

The resulting controller incorporates proportional plus 
integral feedback of the error. This fact is reasonable 
since the system to be controlled, Equation (3.21), is of 
type 0 with a steady state disturbance input. 

To summarize, when LQR synthesis is applied to the 
tracking problem, the proper choice of state variables will 
help to eliminate sensitivity problems and ensure system 
robustness. The state variables should not contain any 
input dependent terms, and they should asymptotically 
approach zero. The system tracking error and its time 
derivative are natural state variable for the LQR synthesis 
procedure. In essence, the key to a successful design 
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depends upon the proper formulation of the performance 
measure. The performance measure of Equation (3.25) 
represents one possible formulation that yields reliable 
results. This performance measure is used in the design of 
the Archytas control systems. These ideas are extended to a 
general multiple input, multiple output (MIMO) system by 
Burl [Ref. 13]. 

D. SUMMARY 

In this chapter, the application of the linear quadratic 
regulator to the tracking problem has been addressed. A 
procedure for formulating a performance measure applicable 
to tracking problems was developed. Additionally, the state 
space representations for both continuous-time and 
discrete-time systems were presented. In the next chapter, 
the results of this chapter are applied to generate control 
systems for Archytas. 
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IV. CONTROL SYSTEM DESIGN FOR ARCHYTAS 



The purpose of this chapter is to use optimal control 
theory to design an automatic flight control system for 
Archytas during the hover mode of flight. Because of the 
coupled nature of the Archytas control problem, a linear 
quadratic regulator (LQR) synthesis technique is pursued. 

A. TOICHYTAS CONTROL SUBSYSTEMS 

The automatic control system is logically separated into 
three subsystems according to the linearized equations of 
Table 2.5. The three control subsystems are: 

1. Roll rate controller; 

2. Altitude rate controller; 

3. Pitch angle and yaw angle controller. 

Because each of these control subsystems is designed 
independently, any cross-coupling which may occur between 
the subsystems is not taken into account. 

B. ARCHYTAS ROLL RATE CONTROLLER 
1. The Roll System 

The roll rate controller for Archytas will serve 
three primary functions: 
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1. Allow the operator to command a desired rotation 
velocity about the vehicle's longitudinal, or x, axis. 
This rotation velocity will permit the operator to 
position Archytas by terminating the rolling motion 
when a desired angle is achieved. This capability is 
critical during the landing phase in order to position 
Archytas correctly with respect to the wind. 

2. Eliminate the rotation velocity imparted to Archytas 
from the effect of cross-winds. 

3 . Eliminate the rotation velocity imparted to Archytas 
from the reactive torques applied to the roll axis as 
the engine speed is varied. 



The simplified equation of motion which describes 
Archytas' roll rate subsystem is given as: 

p=-21.296^ . (4.1) 

The aileron servo dynamics are modeled in Appendix B as a 
second order dynamical system with a natural frequency, , 

of 52.4 rad/sec (8.34 Hz) and a damping ratio, C/ of 0.707. 
The corresponding differential equation is: 

6 =-74 . i6,-2745 . 86 +2745 . 8u„ . (4.2) 

a a a a 

From Equations (4.1) and (4.2), the state equation 
can be written in the matrix form x=Ax+Bu: 



The 



p 




’o 


-21.29 


0 


■p' 




0 




= 


0 


0 


1 






0 


Sa. 




0 


-2745 . 8 


-71.1. 


K. 




2745 . 8. 



roll system tracking error is defined as: 



(4.3) 
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Ep=Pc~P ! 



( 4 . 4 ) 



where is the input command and p is the measured roll 
rate. From Equation (4.4), the differential equation for 
the tracking error is: 

. <<- 5 > 

If Equation (4.5) is combined with the time derivatives of 
Equations (4.1) and (4.2) a new state equation can be formed 
that is appropriate for tracking system design (as discussed 
in Chapter III) . This state equation is: 
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2. Roll Rate Controller Design 

a. Sampling Frequency Selection 

The application of optimal control theory to 
discrete-time systems requires that an appropriate sampling 
frequency be determined. The sampling frequency, f,, is 
simply the inverse of the sampling period, T. A general 
rule of thumb in control systems is to sample a system such 
that 
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( 4 . 7 ) 



where is the maximum bandwidth of the system. From 

Appendix B, the natural frequency of the servos is 52.4 
rad/sec (8.34 Hz). Bassett [Ref. 3: p. 76] demonstrates how 
the Archytas MIMO system can be considered as a set of 
several SISO systems. Each system has a different bandwidth 
for the purpose of determining a sampling frequency. The 
subsequent highest natural frequency is equal to 4.64 
rad/sec (0.74 Hz). Because the natural frequency of the 
servos is a factor of ten greater than the highest natural 
frequency of 4.64 rad/sec, the servo natural frequency will 
be used to determine the required sampling period, T. From 
Equation (4.7), f,=10(52.4 rad/sec) =524 . 0 rad/sec (83.4 Hz), 
using 83.4 Hz as the sampling rate yields 

T=-i- = — 3^=0.012 sec . ( 4 . 8 ) 

fs 83.4 

For the controller designs of this thesis, a sampling period 
of 0.01 seconds will be used. 

b. Discretizing the Roll Rate System 

MATLAB provides tools for computing the discrete- 
time matrix state equation. With the sampling period T 
equal to 0.01 seconds, the discrete-time state space can be 
written in the matrix form x(k+l) x (k) +E U (k) : 



55 



■£p(ic+l)' 
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0 . 001 
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■E^ik)] 
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p(^+l) 




0 1 


-0.2048 


-0 . 0008 


p{k) 




-0 . 0081 


6,(k+i) 
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0 


0 . 8935 


0 . 0067 


5a (^) 
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0.1065 


5jic+i) 




0 


0 


-18 . 5258 


0.3936 . 


5a (^). 




18 . 5253. 



U, {k) . (4.9) 



c. Gain Determination 

The optimal control is determined from the state 
equation and the performance measure: 



00 

k=0 


[Ep{k) p{k) 


6^{k) ^^{k)]Q 


'Ep{k)- 

pik) 

5a (^) 


+ u'^ik) E u ik) ” 








5a {^). 





The optimal control u{k) is: 



u{k) = - 



[^1 ^2 ^3 ^ 4 ] 



\Ep{k) 

pik) 

^,ik) 

^aik) 



= - Ks.{k) 



(4.10) 



(4.11) 



where g is the steady-state gain matrix. The actual control 
u(k) is then obtained by integrating ii{k) to obtain: 



u{k) u(m) = icj ^3 /C 4 ] 



p{m) 

5a (^) 

-k^J2 -k^j){m) -K^h^{m) -k^6^{m) 

m-0 



m=Q 



. 7)=0 



(4.12) 
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The LQR weighting matrices, Q and R, are chosen to satisfy 
the following design criteria: 

1. The overshoot to a step input should be less than five 
percent. 

2. The five percent settling time, tj^, is less than or 
eqpial to 1 second. 



Using an iterative process, the weighting matrices that meet 
the above design criteria were found to be: 



5 = 



0.3 0 0 0 
0 0 0 0 
0 0 0 0 ' 
0 0 0 0 . 



5 =[ 1 ] 



( 4 . 13 ) 



The resulting steady-state feedback gain matrix, K, is: 



= [0.5347 -0.2385 0.1329 0.0017] . ( 4 . 14 ) 

Figure 4.1 shows the roll rate controller block diagram, 
d. Simulation Results 

Figure 4.2(a) shows the response of the closed 
loop system to a step input of six degrees/second. The 
design criteria of an overshoot less than five percent and 
the five percent settling time, tj^, less than one second 
are achieved. Figure 4.2(b) shows the aileron vane angle as 
a result of the step input. Weir [Ref. 16] demonstrates 
that the control vanes stall when displaced by an angle of 
plus or minus 30 degrees from the air flow zero reference. 
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Figure 4.1 Full order Roll Rate 
Controller BlocJc Diagram 
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Time (secs) 

(b) Vane Angle 




Therefore, it is important to ensure that the design not 
require vane angles greater than 30 degrees. 

Because the roll rate controller was designed 
using the linear model, it is necessary to measure the 
relative robustness of the design to ensure that it will be 
applicable to the nonlinear model. The phase and gain 
margins are such measures. From Figure 4.3, the gain margin 
is equal to 32.63 dB and the phase margin is equal to 63.99 
degrees. A general rule of thumb is to require a gain 
margin greater than six dB and a phase margin greater than 
30 degrees to ensure satisfactory performance. Clearly, 
these benchmark values are exceeded and the above design 
should perform well when applied to the nonlinear model. 

The steady-state gains of Equation (4.14) are an 
optimal solution for the performance measure of Equation 
(4.10). However, this optimal solution requires that all of 
the states be measured or estimated. Although measurement 
of the vane angle, 6^, is possible, a computational observer 

would have to be included to provide an estimate of the vane 
velocity, 6^. Inclusion of an observer would add an 

additional task to the onboard digital computer, and 
increase the complexity of the controller design. An 
alternative to the observer would be simply to set the vane 
velocity gain to zero. The system would have to be 
simulated with this gain equal to zero and the phase and 
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Full Order Discrete Open Loop Bode Plot 



Figure 4.3 Full order controller 
Gain and Phase Margins 
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Frequency (rad/sec) 




gain margins computed to determine how the system response 
would be affected. 

e. Reduced Order Model 

Because the servo dynamics are significantly 
faster than all other system dynamics, a second alternative 
is to ignore the servo dynamics completely and form a 
reduced order state space representation of the system. 
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The discrete-time state equation is given as 



'Ep{k^lY 
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( 4 . 16 ) 



Now, the steady-state gain matrix, K/ is determined from 
Equation (4.16) and the performance measure: 



J = 5; j[£p(ic) p(ic) 



ic=0 



Epik) 

.Pik) . 



+ u^(k)E u{k) 



The optimal control ii(k) is: 



u{k) 



~[k, k,] 



\iky 

.Pik) . 



-K p{k) 



( 4 . 17 ) 



( 4 . 10 ) 



and u(k) is 
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(4.19) 



k k 

u(k) = ^ u{im = -V [k.^ k^] 

.73=0 .73=0 

k 



'E^ {m) ' 

i' 

pirn] . 



-k,'^ E^{m) - k^pim) 

.73=0 

The weighting matrices, 2 R, are defined as 



< 2 = 



0.3 0 
0 0 



S=[l] 



( 4 . 20 ) 



The resulting steady-state feedback gain matrix, K, is 

^=[0.4376 -0.2027] . ( 4 . 21 ) 

Figure 4.4(a) shows the response of the system to 
a step input of six degrees/second. The design criteria for 
overshoot and settling time are achieved. The step response 
for this reduced order model is almost identical to the 
response for the full order model given in Figure 4.2(a). 
From Figure 4.4(b), the maximum magnitude of the aileron 
vane angle displacement is 0.47 degrees. This vane 
displacement angle meets the requirement of less than or 
equal to 30 degrees. 

The necessity of determining the relative 
robustness of the system is now magnified due to the reduced 
order controller design. From Figure 4.5, the gain and 
phase margins are equal to 20.59 dB and 57.66 degrees 
respectively. Note that these margins were computed for the 
controller with the full order model. The reduction in gain 
and phase margins from the full order model is an expected 
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Reduced Order Discrete Open Loop Bode Plot 



Figure 4.5 Reduced Order Controller 
Phase and Gain Margins 
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consequence of the reduced order controller. However, the 
requirement for a minimum of six dB gain margin and 20 
degrees of phase margin is satisfied. Therefore, it is 
expected that the reduced order controller design will be 
applicable to the nonlinear model. Figure 4.6 is a 
graphical realization of the reduced order controller 
applied to the roll rate system. 
f . Summary 

An optimal tracker was designed for the roll rate 
control by applying the procedures of Chapter III. This 
control system yielded very satisfactory system performance 
and robustness. Recognizing the dynamics of the servos 
could possibly be ignored, a second reduced order design was 
developed. The resulting controller proved to be robust and 
almost identical in performance to the full order controller 
design. 

Because of the favorable results obtained with 
the reduced order roll rate controller design, the altitude 
rate and the pitch and yaw angle controllers will be 
designed using this technique. The resulting altitude rate 
and pitch and yaw controllers will be evaluated analogously 
to the roll rate controller to ensure they meet or exceed 
the design criteria. 
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Figure 4.6 Reduced Order Roll Rate 
Controller Block Diagram 
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c. 



ARCHYTAS ALTITUDE RATE CONTROLLER 



1. The Altitude System 

The primary function of the altitude rate controller 
for Archytas is to allow uhe operator to command a desired 
translational velocity along the vehicle's longitudinal, or 
X axis. This translational velocity will permit the 
operator to position Archytas at a given altitude by 
terminating the velocity when a desired altitude is 
achieved. This capability is critical during the landing 
phase in order zo control the altitude as Archytas lands in 
a vertical position. 

The simplified equation of motion which describes 
Archytas' altitude rate subsystem is given as: 



Note the change of notation from ii to h , which is more 
appropriate. The throttle servo model is identical to the 
aileron servo model and is given as; 



u = A = 0.10029 6, 



(4.22) 



6j.= -74 . i 6(.-2745 . S6(. + 2745 . 8Ut • 

The Archytas engine is modeled as a first order lag: 



(4.23) 




(4.24) 



The altitude rate tracking error is defined as: 
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( 4 . 25 ) 



,=n„-n ; 



where is the input command and I 2 is rhe measured 

altitude rate. From Equation (4.25), the differential 
equation for the tracking error is: 



E^=-h 



( 4 . 26 ) 



Combining Equation (4.26) with the time derivatives of 
Equations (4.22), (4.23) and (4.24), a state equation that 

is appropriate for tracking system design is: 
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In a manner identical to that used for the roll rate 
controller design, the reduced order altitude rate model is 
formed. The controller is designed using the reduced order 
model and its performance and robustness are evaluated with 
respect to the full order model. Ignoring the servo 
dynamics, the reduced order state space model is: 
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( 4 . 28 ) 
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2. Altitude Rate Controller Design 

a. Discretizing the Altitude Rate System 

With the sampling period, T, equal to 0.01 
seconds, the discrete-time state equation is: 



’ Ef,{k+i) ' 
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b. Gain Determination 

The LQR weighting matrices, Q and R, for the 
altitude rate conrrol loop are chosen to sarisfy ~he 
following design criteria; 

1. The overshoot to a step input should be less than five 
percent. 

2. The five percent settling time, tj^, is less than or 
equal to 2 seconds. 



The resulting weighting matrices are; 



Q= 



0.1 0 0 

0 0 . 01 0 

0 0 0 



[ 1 ] 



(4.30) 



The steady-state feedback gain matrix, K, is: 

^=[-0.30600.20030.0038] . (4.31) 

Figure 4.7 shows the altitude rate controller block diagram. 
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Figure 4.7 Reduced Order Altitude 
Rate Block Diagram 
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c. simulation Results 



Figure 4.3(a) shows the response of the altitude 
rate closed loop system to a step input of ten feet/second. 
The design criteria of an overshoot less than five percent 
and the five percent settling time, less than tv;o 

seconds is achieved. Figure 4.8(b) shows the control 
applied to the system. 

The gain and phase margins were determined as a 
measure of the system robustness. From Figure 4.9, the gain 
margin is equal to 18.03 dB and the phase margin is equal to 
58.78 degrees. These values of gain and phase margin 
indicate that the altitude rate control system design based 
on the reduced order model should provide good results when 
applied to the nonlinear system model. 

D. ARCHYTAS PITCH AND YAW ANGLE CONTROLLER 
1. The Pitch and Yaw Angle system 

The purpose of the pitch and yaw angle controller 
for Archytas is to maintain commanded orientation around the 
pitch and yaw axes. This requirement is necessary to allow 
the operator to position Archytas during landing oy pitching 
or yawing the vehicle slightly to induce a translational 
velocity along the body fixed y or z axis. This control 
problem is complicated by the gyroscopic coupling between 
the axes caused by the propeller. 
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Altitude Rate 



Figure 4.8 Reduced Order Altitude Rate 
Controller step Response 
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Figure 4.9 Reduced Order Altitude Rate 
Controller Gain and Phase Margins 
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The simplified equations of motion that describe the 
pitch and yaw angle subsystem are given as: 



1955-1.62r; 


( 4 . 32 ) 


f = -12 . 286^+1 .78g . 


( 4 . 33 ) 



The elevator and rudder servos are modeled in a manner 
identical to the aileron servo and the corresponding 
differential equations are: 



6^=-74 . 16^-2745 . 86^+2745 . 8U^ ; 

and 


( 4 . 34 ) 


6^=-74 . 16^-2745 . 86^+2745 . 8U^ . 


( 4 . 35 ) 



The pitch and yaw angle tracking errors are defined 

as : 



and 


( 4 . 36 ) 




( 4 . 37 ) 


where and are the input commands and 0 and 7 


are the 



measured pitch and yaw angles, respectively. From Equations 
(4.36) and (4.37), the differential equations for 
the tracking errors are: 



; 


( 4 . 38 ) 


and 
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( 4 . 39 ) 



IT = - r- 
- 

ComJDining Equations (4.38) and (4.39) v/ith the time 
derivatives of Equations (4.34) and (4.35), a state equation 
that is appropriate for tracking system design is x+M Y, 

where the state is: 



x'^ ^ Yq q Eijr Yxn i 0.. J 



( 4 . 40 ) 



the control matrix is ”"=[Ug ii-]", and the system matrices 
are : 
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( 4 . 41 ) 



and 






00000000 2745 . 8 0 

00000000 0 2745 



T 



( 4 . 42 ) 



In a manner similar to that used in the roil rate 
and altitude rate controller designs, the reduced order 
model is formed by ignoring the servo dynamics. The 
resulting state is: 
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( 4 . 43 ) 









jhxa ^ ] 



the control matrix is v~= [6^ 0.]", and the system matrices 
are: 



and 



A = 



0 1 
0 0 
0 0 
0 0 
0 0 
0 0 



0 0 0 0 

-10 0 0 
000 -1.62 
0 10 0 

0 0-1 0 

1.78 0 0 0 



(4.44) 



R 



0 0 -11.19 00 0 

00 0 00 -12.28 



ir 



(4.45) 



2. Pitch and Yaw Angle Controller Design 

a. Discretizing the Pitch and Yaw Angle System 
Using the sampling period, T, equal to 0.01 
seconds, the discrete-time state is: 



x(ic) ^ = [£’e(^) E^ik) q{k) E,,{k) E^{k) (4.46) 

the control matrix is Y.{k) - [6^(A’) and the 



discrete-time system matrices are: 
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(4.47) 



and 



T{k) = 



0 . 0006 
- 0 . 1119 
0 
0 

[ - 0.001 



0 

0 

0 . 001 
0 
0 

- 0 . 1228 



(4.48) 



b. Gain Determination 

The LQR weighting matrices, Q and R, for the 
pitch and yaw angle control loop are chosen to satisfy the 
following design criteria: 

1. The overshoot to a step response should be less than 50 
percent . 

2. The five percent settling time, tj.^, is less than or 
equal to 2 seconds. 



The allowable overshoot requirement may seem too liberal. 
However, due to the small pitch or yaw angles (one to five 
degrees) required to induce translational velocities along 
the body fixed y and z axes, the overshoot should not 
present any problems to the operator. Computer simulations 
have indicated that excessive control values, and vane angle 
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deflections would be required, if the controller was 
designed to limit the overshoot to values less rhan 15 
percent. The large overshoot values are due to the reduced 
order controller. If a full order controller is 
implemented, the large overshoot values could be eliminated 
with small vane angles and a reasonable amount of control. 
However, that would require additional sensors and/or an 
observer as discussed previously. 

The resulting state weighting matrices were found 

to be : 



Q = 



6 0 
0 1 
0 0 
0 0 
0 0 
0 0 



0 0 0 0 
0 0 0 0 
05 0 0 0 

0 6 0 0 
0 0 10 
0 00.05 




The steady-state feedback gain matrix, K, is: 



( 4 . 49 ) 



2.2741 1.8933 -0.6059 0.6390 0.5320 0.0024 ( 4 . 50 ) 

= [-0,6390 -0.5320 -0.0024 2.2741 1.8933 -0.6059J’ 

Figure 4.10 shows the MIMO pitch angle and yaw angle 
controller block diagram [Ref. 14: p. 179]. 
c. Simulation Results 

Figure 4.11(a) show the response of the system to 
a commanded pitch angle of ten degrees and a commanded yaw 
angle of zero. The large value of overshoot was expected 
due to the severe coupling of the axes. Figure 4.11(b) 
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Figure 4.10 Reduced order Pitch Angie and 
Yaw Angle Controller Block Diagram 
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Pitch Angle - Theta 



Figure 4.11 Pitch Angle Equal to Ten 
Degrees / Yaw Angle Equal to Zero 
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Time (secs) 




shows the induced yaw angle from the commanded pitch angle. 
Figure 4.12 gives a similar representation. The yaw angle 
is commanded to five degrees while the pitch angle is held 
at zero. Finally, Figure 4.13 shows the response when the 
pitch and yaw angle are commanded simultaneously to five 
degrees. The maximum vane deflection applied in the 
simulations documented by Figures 4.11, 4.12, and 4.13 was 
1.7 degrees. The design goals for overshoot and settling 
time are satisfied for the linear model with reasonable vane 
deflections. 

d. Singular Value Analysis 

Evaluating the robustness of the roll and 
altitude rate controllers was accomplished by computing 
their respective gain and phase margins. In the case of the 
MIMO pitch and yaw control system, the description of 
robustness in terms of the gain and phase margins becomes 
more complex. Therefore, a different approach is taken in 
order to determine the robustness of the pitch and yaw angle 
control system. The method of choice is to apply singular 
value analysis to the system with perturbations. 

Maciejowski [Ref. 19] presents a detailed development of the 
theory and procedures for computing the singular values of a 
MIMO system and using them for robustness analysis. 
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singular value analysis is used to examine the 
robustness of the MIMO pitch and yaw angle controller. Burl 
[Ref. 13] demonstrates that a system is stable for all input 
multiplicative perturbations, a, such that 



IIaIL i 



1 

Max O X( JW) 
0) 



( 4 . 51 ) 



Where o is the largest singular value of the matrix, and 
T(jcj) is the transfer function from the output of the 
perturbation to the input to the perturbation as shown in 
Figure 4.14. The singular values for the MIMO pitch and yaw 
angle system are plotted in Figure 4.15. The largest 
singular value, is equal to 1.5. Therefore, applying 

Equation (4.51), the MIMO system is expected to be robust 
for perturbations such that the infinity norm of a is less 
than or equal to 0.667. This indicates the controller 
design should provide good results when applied to the 
nonlinear model. [Ref. 13] 



E. RESULTS WITH THE NONLINEAR SYSTEM 

The controllers designed in the previous sections were 
designed using optimal control theory. Because the servo 
states were ignored in determining the gains, the steady- 
state gains are sub-optimal in respect to the full order 
system models. Through simulations, and performance and 
robustness evaluations it was determined that ignoring the 
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Figur* 4.14 MIMO Block Diagram with Perturbations 




o 

< 3 > 



86 




SSV - Pilch/Yaw Angle Control 



Figure 4.15 MIMO Pitch Angle and 
Yaw Angle Singular Values 
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servo dynamics would have no adverse affects. This section 
will apply the controllers designed with respect to the 
linear system models to the nonlinear model. The MATLAB 
programs used to simulate the nonlinear system are included 
as Appendix A. 

1. Simulation One - Figure 4.16 

The system was commanded to climb at a rate of ten 
feet/second for five seconds. The slope of the altitude 
plot indicates that the rate of ten feet/second is achieved. 
The roll rate, pitch angle and yaw angle were set equal to 
zero. The control system is tracking the commanded altitude 
rate input while regulating the other inputs to zero. 

Figure 4.16 shows that the vehicle climbs to an altitude of 
50 feet. Note the coupling between the roll axis (both roll 
rate and roll angle) as the engine accelerates. The 
displacement of the aileron control vanes do not exceed plus 
or minus one degree in deflection. 

2. Simulation Two - Figure 4.17 

Figure 4 . 17 shows the response to a commanded 
altitude rate of 25 feet/second for four seconds. A 
commanded roll rate of ten degrees/second for two seconds is 
applied starting at time equal to six seconds. The slope of 
the altitude plot indicates that the altitude rate of 25 
feet/second is achieved. The vehicle overshoots the desired 
100 feet altitude and begins to approach a steady-state 
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Figure 4.16 Nonlinear Simulation One 
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Figure 4.17 Nonlinear Simulation Two 








90 



value of approximately 80 feet. The large overshoot for the 
altitude is due to the momentum of the vehicle. Climbing at 
a high rate of climb, the vehicle is expected to exhibit 
this type of behavior. Additionally, the control system is 
controlling the altitude rate, not the altitude. The 
operator would have to account for this in his control 
adjustments during flight. Again, the coupling between the 
roll axis and the engine acceleration is observed. The 
commanded roll rate drives the vehicle to an roll angle of 
18°. The two degree error is due to the minus two degrees 
of roll angle caused by the reactive torque to the altitude 
rate. The increased aileron vane displacement in order to 
perform the desired maneuver is shown. 

3. Simulation Three - Figure 4.18 

The vehicle was commanded to climb at a rate of 50 
feet/second for four seconds. A pitch angle of minus five 
degrees was commanded at time equal to six seconds. Figure 
4.18 illustrates that the increased rate of climb increases 
the altitude overshoot. The coupling between the engine and 
roll axis is present. Additionally, the coupling between 
the pitch and yaw axis is shown. As the vehicle is pitched, 
the yaw angle is perturbed. Note, the controller is 
limiting the coupling between the pitch and yaw axis. The 
yaw angle is perturbed less than one degree for a commanded 
pitch angle of minus five degrees. 
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Figure 4.18 Nonlinear simulation Three 





Time (sec) 





92 




4 . Simulation Four - Figure 4.19 

Figure 4 . 19 shows the system response for an 
altitude rate of 50 feet/second for four seconds. The 
commanded roll rate is ten degrees/ second for one second at 
time equal to three seconds. The commanded pitch angle is 
minus three degrees applied at time equal to six seconds. 

The control vanes displacement angles all remain less than 
four degrees for this series of maneuvers. The controller 
is able to provide satisfactory results for the nonlinear 
model . 

F. CONCLUSION 

The controller has demonstrated the ability to control 
the nonlinear system with varying inputs applied. It can be 
concluded that the robustness of each control system is 
sufficient in order to account for the difference in system 
parameters between the linear models and the nonlinear model 
of Archytas. The fact that the control systems maintain 
stability while directing commanded inputs to steer and 
maneuver the vehicle supports the linear modeling approach 
taken in this thesis. The next step in the design process 
is to test the control systems on a prototype vehicle in a 
controlled experiment. Such a test was performed with the 
roll rate control system and an Archytas prototype mounted 
on a test stand. The results of this test will be discussed 
in Chapter V. 
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Figure 4.19 Nonlinear Simulation Four 
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V. CONCLUSIONS 



A. ROLL RATE CONTROL SYSTEM FIELD TEST 

Figure 5.1 shows the Archytas prototype (minus the wings 
and canard) mounted on a test stand that allows the vehicle 
to spin about the longitudinal, or x axis. This 
configuration was used during the testing of the roll rate 
controller. The testing was accomplished in a qualitative 
manner; no empirical data was collected, the object of the 
test being to validate the roll rate controller design. The 
test consisted of three parts: 

1. The ability of the roll rate controller to eliminate 
the rotational velocity imparted to Archytas from the 
reactive torques applied to the roll axis as the engine 
speed is varied (decoupling) . 

2. The ability of the roll rate controller to eliminate 
the rotational velocity imparted to Archytas from the 
effect of cross-winds (disturbance rejection) . 

3. The ability of the roll rate controller to allow the 
operator to impart or terminate a rotational velocity 
about the x axis (tracking commanded inputs) . 



Part one was accomplished by varying the engine speed 
from idle to maximum rpm with the commanded roll rate equal 
to zero. Each time the engine speed was varied, the 
rotational velocity of the vehicle was eliminated by the 
roll rate controller. It was observed that the rotational 
velocity was eliminated within approximately two seconds. 
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Figure 5.1 Archytas Prototype Mounted 
on the Test Stand 
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This response is consistent with both the linear and 
nonlinear simulations. 

Part two was accomplished by setting the engine speed at 
constant values ranging from idle to maximum rpm with the 
commanded roll rate equal to zero. A disturbance about the 
X axis was imparted to Archytas by pushing the test stand 
mounting bracket connected to the vehicle. The effect of 
pushing the connecting bracket is analogous to a rotational 
velocity imparted by cross-winds. It was observed that the 
rotational velocity was eliminated within approximately two 
seconds. This response is consistent with both the linear 
and nonlinear simulations. 

Part three was accomplished by inputting a commanded 
roll rate while the engine speed was held constant and also 
varied. It was observed that the operator could position 
Archytas at a desired angle under either test condition. 
Additionally, disturbances were imparted to Archytas during 
these tests. The roll rate controller proved robust enough 
to eliminate the disturbances and allow the operator to 
position the vehicle. This response is consistent with both 
the linear and nonlinear simulations. 

Based on the computer simulations, both linear and 
nonlinear, and the results of the above tests, it is 
concluded that the roll rate controller design is valid. 
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B. FUTURE RESEARCH 



The development of the Archytas control systems was 
accomplished using the data from the AROD as first 
approximations. Future research will include the 
computation of Archytas specific data through empirical 
methods or obtaining the data by wind tunnel testing. This 
would allow for a better model of the Archytas to be 
developed. 

The second order dynamical model of the servos was 
obtained with the vanes mounted on the servos with zero 
downwash from the propeller. Future research should include 
modeling of the servos with the engine at hover rpm and the 
vanes positioned within the downwash. 

The effects of inclusion of the vane position angles in 
the control systems could be evaluated further. In 
particular, the effect the vane angle position feedback 
would have in the reduction of the overshoot in the pitch 
and yaw angle controller should be investigated. 
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APPENDIX A 



MATLAB SIMULATION PROGRAMS 






‘S'o'o'o'b'o'o'o'b'b'b'b'b'b'o'b'b'o'b'b'b'o'o'b'o'o 






:%%s 



% PROGRAM NAME: nonlin.m 

% 

% DESCRIPTION: 

% This program will simulate the Archytas nonlinear model 
% with the linear controllers. The Control Laws are the 
% Steady-state Linear Quadratic Regulator solutions with 
% the necessary modifications to the performance measure. 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S’S-S-S-S-S-S-S-?-?-?-?-?-?-?-?-?-?-?-?-?- 

•S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'o'o'&'o'S'S'S'&'S’S’S'S'O’S'o'S'^'o'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S 

clear 

% Call the program noninp.m to enter the desired roll rate, 

% alttitude rate, pitch and yaw angles. 

noninp 

% Call the program in_cond.m to load the initial conditions 

and % constants for the simulation. 

in_cond 

% Call the program n_l_inp.m to load the wind tunnel data. 
n_l_inp 

% Loop for simulation 
for k=l : num_of_stepsl 



% Compute the speed of the propeller 
speed (k) = speed_hover + del_rpm(k) ; 

% Compute the thrust supplied by the engine to the vehicle 
thrust(k) = thrust_hover + Xrpm * del_rpm(k) ; 

% Check thrust limits 
if thrust (k) < thrust_min 
thrust (k) = thrust_min; 
elseif thrust (k) > thrust_max 
thrust (k) = thrust_max; 

end 



velocity_tip (k) = speed_hover + del_rpm(k) ; 
del_tip(k) = velocity _tip (k) - speed_hover ; 

tip_del(k) = velocity_tip (k) - speed (k) ; 

stad(k) = tip_del(k); 

% Compute vehicle velocity and angle of attack total 
aoatot(k) = alpha_max; 

velocity_tot(k)=sqrt(u(k) ^2+v(k) ^2+w(k) '' 2 ) ; 
if velocity_tot (k) > vaero 

vwterm = sqrt(v(k)''2 + w(k)^2); 
aoatot(k) = asin (vwterm/velocity_tot (k) ) ; 
if aoatot(k) < alpha_min 
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o\o o\o o\o o\o o\o o\o o\o o\o 



aoatot(k) = alpha_min; 
elseif aoatot(k) > alpha_max 
aoatot(k) = alpha_max; 

end 

end 



veq(k) = tablel (veq_table, aoatot (k) ) ; 
newveq(k) = sqrt (weight_ratio) *veq (k) ; 
rhova(k) = rho * newveq(k) * pi; 
velocity _tip_2 (k) = 2 . 0/velocity_tip (k) ; 
velocity_delta (k) = velocity_tot (k) - newveq(k); 



% This section assumes that no aerodynamic forces and moments 
% exist 

if velocity_tot (k) <= vaero 



fax 

fay 

faz 

mrx 

mpy 

myz 

pitch_trim 

yaw_trim 

pitch_factr 

yaw_factr 



= thrust (k) 

= 0 . 0 ; 

= 0 . 0 ; 

= 0 . 0 ; 

= 0 . 0 ; 

= 0 . 0 ; 

= 0 . 0 ; 

= 0 . 0 ; 

= 1 . 0 ; 

= 1 . 0 ; 



* gravity /weight; 



elseif velocity _tot (k) > vaero 
if aoatot (k) <= alpha_min 



fax 


=r 


thrust (k) 


fay 


= 


0.0; 


faz 


= 


0.0; 


mrx 


= 


0.0; 


mpy 


= 


0.0; 


myz 


= 


0.0; 


pitch_trim 


= 


0.0; 


yaw_trim 


= 


0.0; 


pitch_factr 


= 


1.0; 


yaw factr 


= 


1.0; 



end 



* gravity /weight; 



else 



% This section computes the aerodynamic forces and moments 

lleq = weight_ratio * leg; 

ddeq = weight_ratio * deq; 

sseq = weight_ratio * seqq; 

req = tablel (req_table, aoatot (k) ) ; 
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peq 

yeq 



= tablel (peq_table, aoatot (k) ) ; 
= tablel (yeq_table, aoatot (k) ) 

rreq = weight_ratio- * req; 
ppeq = weigth_ratio * peq; 
yyeq = weight_ratio * yeq; 



cldel = tablel (cldel_table, aoatot (k) ) ; 

cddel = tablel (cddel_table, aoatot (kn / 

csdel = tablel (csdel_table, aoatot (k) ) ; 

Idel = velocity_tip_2 (k) *lleq-rhova*weight_ratio*cldel ; 
ddel = velocity_tip_2 (k) *ddeq-rhova*weight_ratio*cddel; 
sdel = velocity_tip_2 (k) *sseq-rhova*weight_ratio*csdel ; 



crdel = tablel (crdel_table, aoatot (k) ) ; 
cpdel = tablel (cpdel_table, aoatot (k) ) ; 

cydel = tablel (cydel_table, aoatot (kn ; 

rdel = velocity_tip_2 (k) *rreq-rhova*weight_ratio*crdel ; 
pdel = velocity_tip_2 (k) *ppeq-rhova*weight_ratio*cpdel ; 
ydel = velocity_tip_2 (k) *yyeq-rhova*weight_ratio*cydel ; 



Islope 

dslope 

sslope 

ffl 

ffd 

ffs 



rslope 

pslope 

yslope 



tablel ( lslope_ 
tablel (dslope 
tablel (sslope" 



table, aoatot (k) ) ; 
table, aoatot (k) ) ; 
table, aoatot (k) ) ; 



lleq+ldel*del_tip (k) +lslope*weight_ratio* 
velocity_delta (k) ; 

ddeq+ddel*del_tip (k) +dslope*weight_ratio* 
velocity_delta (k) ; 

sseq+sdel*del_tip (k) +sslope*weight_ratio* 
velocity_delta (k) ; 

tablel (rslope_table, aoatot (k) ) ; 
tablel (pslope_table, aoatot (k) ) ; 
tablel (yslope_table, aoatot (k) ) ; 



mr = rreq+rdel*del_tip (k) +rslope*weight_ratio* 

velocity_delta (k) ; 

mp = ppeq+pdel*del_tip (k) +pslope*weight_ratio* 

velocity _delta (k) ; 

my = yyeq+ydel*del_tip (k) +yslope*weight_ratio* 

velocity_delta (k) ; 

fs = ffs * gravity /weight ; 

fl = ffl * gravity /weight; 

fd = ffd * gravity / weight ; 

delta = atan(v(k) /w(k) ) ; 



fax = fl * sin(aoatot (k) ) - fd * cos (aoatot (k) ) ; 
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fayz 


=-fl * 


cos (aoatot (k) ) 


- fd * 


sin (aoatot (k) ) ; 


fay 


= fayz 


* sin (delta); 






f az 


= fayz 


* cos (delta); 






mrx 


= my * 


sin (aoatot (k) ) 


- mr * 


cos (aoatot (k) ) ; 


mpyyz 


=-my * 


cos (aoatot (k) ) 


- mr * 


sin (aoatot (k) ) ; 


mpy 


= mpyyz 


* sin (delta); 






myz 


= mpyyz 


* cos (delta); 







velocity_trim = tablel (velocity_triin_table, aoatot (kl) ) ; 
vaneff = tablel (vaneff table, aoatot (k) ) ; 



pitch_trim = velocity_trim * cos (delta); 

yaw_trim =-velocity_trim * sin(delta); 

pitch_factr = vaneff; 

yaw_factr = vaneff; 



end 

lat(k) = ( (rscale*mrx*gravity*144 . 0) /Ixx) + (La*del_a (k) ) - 
(velocity_delta (k) ^2 /Ixx) * prop_torq; 

mat(k) = ( (mpy*gravity*144 . 0) /lyy) + (Me*pitch_f actr* 

(del_e (k) -pitch_trim) ) ; 

nat(k) = ( (myz*gravity*144 . 0) /Izz) + (Nr*yaw_f actr * 

(del_r (k) -yaw_trim) ) ; 

% Begin Differential Equations of Motion 

% Altitude Differential Equations 
tempi (k) =u (k) *cos (theta (k) ) *cos (psi (k) ) ; 
temp2 (k) =v(k) * (sin (phi (k) ) *sin (theta (k) ) *cos (psi (k) ) - 
cos (phi (k) ) *sin (psi (k) ) ) ; 

temp3 (k)=w(k) *(cos(phi(k) ) *sin(theta(k) ) *cos(psi(k) )+ 
sin (phi (k) ) *sin (psi (k) ) ) ; 
u_earth(k)= tempi (k) +temp2 (k) +temp3 (k) ; 

% Heading Angle Differential Equations 
temp5 (k) =p (k) *cos (theta (k) ) *cos (psi (k) ) ; 
temp6 (k) =q (k) * (sin (phi (k) ) *sin (theta (k) ) *cos (psi (k) ) - 
cos (phi (k) ) *sin (psi (k) ) ) ; 

temp? (k) =r (k) * (cos (phi (k) ) *sin (theta (k) ) *cos (psi (k) ) 

+sin (phi (k) ) *sin (psi (k) ) ) ; 
hdg_dot (k) =temp5 (k) +temp6 (k) +temp7 (k) ; 

% Euler Angle Differential Equations: Pitch, Yaw, and Roll 
theta_dot (k) =q (k) *cos (phi (k) ) -r (k) *sin (phi (k) ) ; 
phi_dot=p (k) + (q (k) *sin (phi (k) ) +r (k) *cos (psi (k) ) ) * 
tan (theta (k) ) ; 
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psi_dot (k) = (q (k) * sin (phi (k) ) +r (k) *cos (phi (k) ) ) / 
cos (theta (k) ) ; 

% Vehicle Pitch, Yaw and Roll Rate Differential Equations WRT 
% Body Axes 

p_dot(k) = (bca*q (k) *r (k) ) - (Iralph*stad (k) ) +lat (k) ; 
q_dot(k) = (cab*p (k) *r (k) ) -(Irx*r (k) *speed (k) /lyy) +mat (k) ; 
r_dot(k) = (abc*p(k) *q(k) )+(Irx*q(k) *speed(k) /Izz) +nat (k) ; 

% Vehicle Velocity Differential Equations WRT Body Axes 
u_dot (k)=( (v(k)*r(k) )-(w(k)*q(k) ) ) -gravity *cos (theta (k) ) * 
cos (psi (k) ) +fax; 

v_dot (k) = ( (w(k) *p (k) ) - (u(k) *r (k) ) ) +gravity*cos (theta (k) ) * 
sin (psi (k) ) +fay ; 

w_dot (k) = ( (u (k) *q (k) ) - (v (k) *p (k) ) ) -gravity* 
sin (theta (k) ) +faz ; 



% Define Error Variable in Pitch and Yaw 
E_theta(k)=theta_com(k) -(theta (k)+noise(k) ) ; 

E_psi (k) =psi_com(k) - (psi (k) +noise (k) ) ; 

% Define Error Variable in Roll 
E_p (k) =p_com (k) - (p (k) tnoise (k) ) ; 

% Define Error Variable in Altitude 
E_h_dot (k) =h_dot_com(k) - (u(k) +noise (k) ) ; 

% Integrate Errors 

E_thetain (k+1) =E_thetain (k) +Ts*E_theta (k) ; 

E_psiin (k+1) =E_psiin (k) +Ts*E_psi (k) ; 
E_pin(k+l)=E_pin(k)+Ts*E_p(k) ; 

E_h_dotin (k+1) =E_h_dotin (k) +Ts*E_h_dot (k) ; 

% Define Control System Gains 
% Roll Rate Control Gains 
Kl=[ 0.4376 -0.2027]; 

% Pitch and Yaw Angle Control Gains 

K2=[ 2.2741 1.8933 -0.6059 0.6390 0.5320 0.0024; 

-0.6390 -0.5320 -0.0024 2.2741 1.8933 -0.6059]; 

% Altitude Rate Control Gains 
K3=[-0.3060 0.2003 0.0038]; 

% Calculate the Aileron (roll) Servo Control Input 
Ua(k)=-Kl*[E_pin(k+l) ; (p (k) +noise (k) ) ] ; 

% Calculate the Elevator (pitch) Servo Control Input 
Ue(k)=-K2(l, : ) * [E_thetain (k+1) ;E_theta(k) ;q(k) ;E_psiin (k+1) ; 

E_psi(k) ;r(k) ] ; 
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% Calculate the Rudder (yaw) Servo Control Input 

Ur (k)=-K2 (2, : ) * [E_thetain(k+1) ;E_theta(k) ;q(k) ;E_psiin(k+l) ; 

E_psi(k) ;r (k) ] / 

% Calculate the Throttle (altitude) Servo Control Input 
Ut(k)=-K3*[E_h_dotin(k+l) ;u(k) ;del_rpiti(k) ] ; 

% Begin Integration of Equations of Motion 
del_a_dot_dot (k)=-Hl*del_a_dot (k) -H2*del_a(k) +H2*Ua(k) ; 
del_e_dot_dot (k) =-Hl*del_e_dot (k) -H2*del_e (k) +H2*Ue (k) ; 
del_r_dot_dot (k) =-Hl*del_r_dot (k) -H2*del_r (k) +H2*Ur (k) ; 
del_t_dot_dot(k)=-Hl*del_t_dot(k) -H2*del_t(k) +H2*Ut(k) ; 

%Roll, pitch, yaw 
p(k+l)=p(k) +Ts*p_dot(k) ; 
q (k+1) =q (k) +Ts*q_dot (k) ; 
r (k+l)=r (k) +Ts*r_dot (k) ; 

%Velocities 

u(k+l) =u(k) +Ts*u_dot (k) ; 

V (k+1) =v (k) +Ts*v_dot (k) ; 
w(k+l)=w(k)+Ts*w_dot(k) ; 

%Euler angles 

theta (k+1) =theta (k) +Ts*theta_dot (k) ; 
phi (k+1) =phi (k) +Ts*phi_dot (k) ; 
psi (k+1) =psi (k) +Ts*psi_dot (k) ; 

%Servos 

del_a_dot (k+1) =del_a_dot (k) +Ts*del_a_dot_dot (k) ; 
if del_a_dot (k+1) > rmax 
del_a_dot (k+1) = rmax; 
elseif del_a_dot (k+1) < -rmax 
del_a_dot (k+1) = -rmax; 

end 

del_e_dot (k+1) =del_e_dot (k) +Ts*del_e_dot_dot (k) ; 
if del_e_dot (k+1) > rmax 
del_e_dot (k+1) = rmax; 
elseif del_e_dot (k+1) < -rmax 
del_e_dot (k+1) = -rmax; 

end 

de l_r_dot ( k+ 1 ) =de l_r_dot ( k ) +Ts *de l_r_dot_dot ( k ) ; 
if del_r_dot (k+1) > rmax 
del_r_dot (k+1) = rmax; 
elseif del_r_dot (k+1) < -rmax 
del_r_dot (k+1) = -rmax; 

end 
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del_t_dot (k+1) =del_t_dot (k) +Ts*del_t_dot_dot (k) ; 
if del_t_dot (k+1) > rmaxt 
del_t_dot (k+1) = rmaxt; • 
elseif del_t_dot (k+1) < -rmaxt 
del_t_dot (k+1) = -rmaxt; 

end 

% Engine 

del_rpm_dot (k+1) =-w_e*del_rpm(k) +Ke*w_e*del_t (k) ; 

del_a (k+1) =del_a (k) +Ts*del_a_dot (k) ; 

del_e (k+1) =del_e (k) +Ts*del_e_dot (k) ; 
if del_e(k+l) > maxdfl 
del_e(k+l) = maxdfl; 
elseif del_e(k+l) < -maxdfl 
del_e(k+l) = -maxdfl; 

end 

del_r (k+1) =del_r (k) +Ts*del_r_dot (k) ; 
if del_r(k+l) > maxdfl 
del_r(k+l) = maxdfl; 
elseif del_r(k+l) < -maxdfl 
del_r(k+l) = -maxdfl; 

end 

del_t (k+1) =del_t (k) +Ts*del_t_dot (k) ; 
if del_t(k+l) > tamax 
del_t(k+l) = tamax; 
elseif del_t(k+l) < -tamax 
del_t(k+l) = -tamax; 

end 

% Engine 

del_rpm(k+l) =del_rpm(k) +Ts*del_rpm_dot (k) ; 

% Altitude 

alt (k+1) =alt (k) +Ts*u_earth (k) ; 

% Heading Angle 

heading (k+1) =heading (k) +Ts*hdg_dot (k) ; 

t(k)=Ts* (k-1) ; 

% End of simulation loop 
end 

t(k+l)=Ts*k; 

u_earth (k+1) =u earth (k) ; 
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% Plot the Results of the Simulation 

subplot (221) ,plot(t,alt) ,grid,xlabel('Time 
(sec) ' ) ,ylabel( 'Altitude (ft) ' ) ; 
title ( 'Altitude' ) ; 

subplot(222) ,plot(t,p*57.2958) , grid, xlabel( 'Time 
(sec) ') ,ylabel ( 'Roll Rate (deg/sec)'); 
title('Roll Rate - p'); 

subplot(223) ,plot(t,phi*57.2958) ,grid,xlabel('Time 
(sec) ') ,ylabel( 'Angle (deg)'); 
title ('Roll Angle Phi'); 

subplot(224) , plot (t, theta*57 . 2958 ) , grid , xlabel (' Time 

(sec) ') ,ylabel ( 'Angle (deg)'); 

title ('Pitch Angle Theta'); 

meta nonl 

subplot ( 111) ; 



subplot(221) ,plot(t,psi*57.2958) ,grid,xlabel('Time 
(sec) ') ,ylabel( 'Angle (deg)'); 
title ('Yaw Angle Psi'); 

subplot (222) ,plot(t,del_a*57.2958) ,grid,xlabel('Time 
(sec) '), ylabel ( 'Angle (deg)'); 
title ( 'Aileron Vane Angle'); 



subplot(223) , plot ( t , de l_r *5 7 . 2 9 5 8 ) , 
(sec) ') ,ylabel( 'Angle (deg)'); 
title ( 'Rudder Vane Angle'); 
subplot(224) , plot ( t , del_e*5 7 . 2 9 5 8 ) , 
(sec) '), ylabel ( 'Angle (deg)'); 
title ( 'Elevator Vane Angle'); 
meta non2 
subplot ( 111) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



grid,xlabel( 'Time 
grid, xlabel ( ' Time 



%%% 



• ^ ^ S' 5^ ^ ^ ^ ^ ^ ^ 






% PROGRAM NAME: noninp.m 



o'o'o'o'o'a'o'o' 






•6 

% 



DESCRIPTION: 

This program prompts the user for the desired 
simulation inputs. 



S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'SrS'S'S'^^^ 












*0 'o *0 'b 






% Prompt the user for the desired simulation time. 
disp(' Enter the simulation time in seconds.'); 
simulation_time = input('>>>> '); 

% Prompt the user to enter to desired sampling time - Ts. 
disp(' Enter the desired sampling time - Ts.'); 
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o\o o\o o\o o\o o\o 



Ts = input ('>>>> '); 

% Compute the number of simulation steps based on the 
desired simulation time. 

num_of_stepsl = round (simulation_time/Ts) ; 

% Prompt the user for the length of the altitude rate step 
input time. 

disp(' Enter the length of the altitude rate input in 
seconds. ' ) ; 

step_input_timel = input ('»>> '); 

% If step_input entered is greater than simulation_time 
entered, set them equal, 
if step_input_timel > simulation_time 
step_input_timel = simulation_time; 

end 

if step_input_timel -= 0 

% Prompt the user for the desired magnitude of the step 
% input in # of deg/ sec. 

disp(' Enter the magnitude of the altitude rate input in 
feet/sec . ' ) 

h_dot = input ( ' »» ' ) ; 
elseif step_input_timel == 0 
h_dot = 0.0; 

end 

% Compute the number of step input steps based on the 
% desired step input time. 
num_of_steps2 = step_input_timel/Ts; 

num_of_steps2 = round (num_of_steps2 ) ; 

% Compute the number of zeros to be padded to the input 
% based on the length simulation time and length of the 
% desired step input time. 

num_of_zerosl = num_of_stepsl - num_of_steps2 ; 

% Define the commanded input vector 
h_dot_com = [h_dot * ones ( 1 , num_of_steps2 ) 
zeros ( 1 , num_of_zerosl) ] ; 

% Prompt the user for whether he wants to include a roll 
rate in the simulation. 

ans = input ('Do you want to include a roll rate in the 
simulation? y/n [y] : ','s'); 
if isempty(ans) 
ans ='y' ; 

end 

if ans =='y' 
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% Prompt the user for the start time for the roll rate 
step input. 

disp(' Enter the start time for the roll rate input.'); 
start_timel = input ('>>>> '); 

% Prompt the user for the length of the roll rate input. 
disp(' Enter the length of the roll rate imput in 
seconds . ' ) ; 

length_roll_rate = input ('>>>> '); 

% Compute the number of zeros to pad front of the 
% command vector. 
num_of_zeros2 = start_timel/Ts; 
num_of_zeros2 = round ( num_of_zeros2 ) ; 

% Computer the number of step input steps and number of 
% zeros to pad back of the command vector 

num_of_steps3 = length_roll_rate/Ts; 
num_of_steps3 = round (num_of_steps3 ) ; 
num_of_zeros3 = num_of_stepsl- 

round( (start_timel+length_roll_rate) /Ts) ; 

% Prompt the user for the desired magnitude of the roll 
% rate input in # deg/se. 

disp(' Enter the magnitude of the roll rate in 
degrees/ sec. ') ; 
pi = input ('>>>> '); 

% Define the commanded input vector 

p_com = [ zeros (l,num_of_zeros2) pi * 0.017453 * 

ones ( 1 , num_of_steps3 ) zeros ( 1 , num_of_zeros3 ) ] ; 

else 

p_com = [ zeros ( 1, num_of_stepsl) ] ; 

end 

% Prompt the user for whether he wants to include a pitch 
% angle in the simulation. 

ans = input ('Do you want to include a pitch angle in the 
simulation? y/n [y] : ','s'); 
if isempty(ans) 
ans = 'y ' ; 

end 

if ans =='y' 

% Prompt the user for the start time for the pitch angle 
% command. 

disp(' Enter the time at which you desire to input the 
pitch angle.') 

start_time2 = input ('>>>> '); 

% Compute the number of zeros to pad front of the 
% command vector 

num_of_zeros4 = round (start_time2/Ts) ; 
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% Compute the number of ones to multiply the desired 
% angle by to create the command vector. 

number_of_onesl = num_of_stepsl - round (start_time2/Ts) ; 
% Prompt the user for the desired pitch angle in 
% degrees. 

disp(' Enter the desired pitch angle in degrees.'); 
theta_com = input ('»» '); 

theta_com = [ zeros (l,num_of_zeros4) theta_com * 0.017453 
* ones ( 1 , number_of_onesl ) ] ; 

else 

theta_com = [ zeros ( 1, num_of_stepsl) ] ; 

end 



% Prompt the user for whether he wants to include a yaw 
% angle in the simulation. 

ans = input ('Do you want to include a yaw angle in the 
simulation? y/n [y] : ','s'); 

if isempty(ans) 
ans = 'y ' ; 

end 

if ans =='y' 

% Prompt the user for the start time for the yaw angle 
% command. 

disp(' Enter the time at which you desire to input the 
yaw angle. ' ) 

start_time2 = input ('>>» '); 

% Compute the number of zeros to pad front of the 
% command vector 

num_of_zeros4 = round ( start_time2 /Ts ) ; 

% Compute the number of ones to multiply the desired 
% angle by to create the command vector. 
number_of_onesl = num_of_stepsl - round (start_time2/Ts) ; 
% Prompt the user for the desired pitch angle in 
% degrees. 

disp(' Enter the desired yaw angle in degrees.'); 
psi_com = input ('>»> '); 

psi_com = [ zeros (l,num_of_zeros4) psi_com * 0.017453 * 
ones ( 1 , number_of_onesl ) ] ; 

else 

psi_com = [ zeros ( 1, num_of_stepsl) ] ; 

end 

% Prompt the user for whether he wants to add sensor noise 
% to the simulation. 

ans = input ('Do you want to include sensor noise in the 
simulation? y/n [y] : ','s'); 
if isempty(ans) 
ans = ' y ' ; 

end 

if ans =='y' 
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% Define the random vector to represent the measurement 
% noise. 

% Assume the accuracy is +/- 1.745e-04. 
rand ( 'uniform' ) ; 
noise = rand (l , num_of_stepsl) 
vard = (1/12) * (1.745e-04) 

A = sqrt (vard/cov (noise) ) 



2 



Desired variance 
Scalar multiplier to 
change the variance 



noise = A .* noise; 

% Adjust the mean to zero 
noise = noise - mean (noise) ; 
else 

noise = [ zeros ( 1 , num_of_stepsl) ] ; 

end 






> o o o 



o o o o o 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' 



PROGRAM NAME: in cond.m 



% 

% DESCRIPTION: 

% This program inputs the initial conditions and 
% constant values used in nonlin.m 



> *6 'S iS *6 



) *S 'b 'o *o *b *0 



'o'o'o'o'b'o'o'o'o'o'ob'^'o'o'o'o'o'b'o'o^'o'o'o'o'o'o'b 



“O ^ ^ O O "O TO O 70 TO 






% CONSTANTS 

vaero=0.5; speed_max=859 . 0 ; speed_min=649 . 0 ; 
prop_torq=0 . 0729 ; thrust_min=35. 0; thrust_max=105 . 0 ; 
tamax=100.0; rmaxt=100.0; 
maxdf 1=0 . 5236 ; rmax=0 . 87266 ; 

La=-21.04; Me=-9.01; Nr=-11.40; 

Hl=71.1; H2=2745.8; 

w_e=2.0; Ke=837.758; Xrpm=0.2122; speed_hover=712 . 0943 ; 
alpha_min=0 . 174533 ; alpha_max=l . 570796 ; 
gravity=32 . 174 ; pi=3 . 1415962 ; rad_to_deg=180 . 0/pi; 
Ixx=6908.4; Iyy=22944 . 64; Izz=21515.64; Irx=41.62; 
bca=0. 20685; cab=0. 63663; abc=-0 . 74533 ; Iralph=0.00916; 
weight=85.0; thrust_hover=85 . 0 ; velocity_tip=722 . 5663 ; 
rho=0. 00192; weight_ratio=l . 0 ; 

leq=85.0; deq=0.0; seqq=0.0; sslope=0 • 0801 ; rscale=0.0; 
% INITIAL CONDITIONS 



p(l)=0.0; 

q(l)=0.0; 

r(l)=0.0; 
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o\o o\o o\o o\o o\o 



u(l)=0.0; 

v(l)=0.0; 

w(l)=0.0; 



phi(l) =0.0; 
psi(l) =0.0; 
theta ( 1) =0 . 0 ; 



alt(l) =0.0; 
dist(l) =0.0; 
heading (1) =0.0; 



0 . 0 ; 
0 . 0 ; 
0 . 0 ; 
0 . 0 ; 

del_a(l) =0.0; 
del_e(l) =0.0; 
del_r(l) =0.0; 
del_t(l) =0.0; 
del_rpm(l)=0.0; 



E_thetain(l) = 
E_psiin(l) = 
E_pin(l) = 
E_h_dotin(l) = 



del_a_dot ( 1) =0 . 0 ; 
del_e_dot ( 1 ) =0 . 0 ; 
del_r_dot ( 1 ) =0 . 0 ; 
del_t_dot ( 1) =0 . 0 ; 



;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 



o o o o c c o 
o o o o o o o 






PROGRAM NAME: 



n 1 inp . m 



-6 



DESCRIPTION: 

This file inputs the wind tunnel data into the 
MATLAB wordspace 



^%%%%%%% 

>0000000 



aoatotl=[0. 174533 0.349066 0 
0.959931 1.047198 1 
1.570796] ; 



CRdell =[0.00 


0.00 


0.03 


0.06 


0.07 


0.08 


0.07 


0.07 


Rslopel=[0. 00 


0.01 


0.07 


0.14 


0.18 


0.18 


0.18 


0.17 



523599 0.785398 0.872665... 
134464 1.221730 1.308997... 



0.06 

0.07 0.05]; 
0.16 

0.18 0.20]; 
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oV> o\o o\o 



Reql =[2.16 4.84 6.66 7.23 7.76 

8.66 8.92 6.89 6.66 5.68 0.00]; 



CPdell =[-0.08 -0.09 -0.08 -0.02 -0.00 

0.03 0.07 0.12 0.08 0.12 0.00]; 

Pslopel=[-0.29 -0.24 -0.21 -0.04 -0.01 

0.08 0.17 0.29 0.19 0.31 0.40]; 

Peql =[-5.06 2.24 7.46 15.32 16.65 17.12 
17.37 17.21 18.03 16.10 0.00]; 

CYdell =[0.03 0.04 0.04 0.03 0.03 

0.03 0.03 0.01 0.01 0.01 0.00]; 

Yslopel=[0. 12 0.12 0.11 0.08 0.07 

0.06 0.06 0.02 0.02 0.02 0.00]; 

Yeql =[4.54 6.36 6.82 3.13 3.05 

3.73 3.17 2.34 2.34 1.81 1.00]; 

VanEf f 1=[ 1. 50 1.45 1.40 1.35 1.30 

1.25 1.20 1.15 1.10 1.05 1.00]; 

Veql =[172.4 115.4 87.65 62.85 55.45 

48.00 41.2234.68 28.98 22.37 0.0]; 

CLdell =[0.16 0.28 0.32 0.29 0.30 0.26 

0.23 0.19 0.20 0.16 0.00]; 

Lslopel=[0.60 0.81 0.82 0.73 0.74 0.64 

0.56 0.46 0.50 0.40 0.00]; 

CDdell =[0.44 0.46 0.46 0.48 0.50 0.50 

0.50 0.50 0.54 0.52 0.50]; 

Dslopel=[l. 67 1.32 1.19 1.20 1.23 1.22 

1.24 1.23 1.34 1.32 1.25]; 

CSdell =[-0.01 -0.02 -0.03 -0.01 -0.01 

-0.02 0.01 0.00 0.02 0.04 0.06]; 

vtriml =[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] 



crdel_table =[aoatotl'CRdell']; 
rslope_table=[aoatotl'Rslopel' ] ; 
req_table =[aoatotl' Reql']; 



cpdel_table =[aoatotl' 
pslope_table= [ aoatotl ' 
peq_table =[ aoatotl' 



CPdell' ] ; 
Pslopel ' ] ; 
Peql' ] ; 



cydel_table =[ aoatotl' 
yslope_table= [ aoatotl ' 
yeq_table =[ aoatotl' 



CYdell' ] ; 
Yslopel' ] ; 
Yeql' ] ; 



vanEff_table=[ aoatotl' VanEf fl' ] ; 
veq_table =[ aoatotl' Veql']; 
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cldel_table =[aoatotl' CLdell']; 
lslope_table= [ aoatot 1 ' Lslopel ' ] 

cddel_table =[ aoatot 1' CDdell']; 
dslope_table=[ aoatot Dslopel'] 

csdel_table =[aoatotl' CSdell']; 
velocity_trim_table =[ aoatot 1' 



f 

Vtriml' ] 



APPENDIX B CONTROL SERVOS 



The control servos incorporated into Archytas are 
position-commanded, quarter-scale model airplane servos. 

Five of these servos are used on Archytas: four to drive the 
aerodynamic control surfaces and one to drive the engine 
throttle. These servos can be modeled as second-order 
dynamical systems; 



where C(s) is the output and R(s) is the input. ^ and 
are referred to as the damping ratio and the natural 
frequency, respectively. 

The values for f and were determined by experiment. 
The control vanes were mounted on the servos and position 
data was obtained for a range of step input commands. The 
input commands ranged from three to thirty degrees of vane 
deflection. By plotting the servo angle response data 
versus time, a series of step input response curves was 
developed. Using MATLAB, trial and error was used to fit a 
second order prototype system to the step response curves of 
the servos. Figure B.l shows the "best fit" model. This 
model is defined by ^=0.707 and cj„=52.4 radians/second. 



C(s) 

R(s) 




(B.l) 






2 

n 
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Figure B.l Servo Response curve 
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Thus, the second order servo model used in the design of the 
Archytas control systems is: 

CM. = 2745.8 

R{s) s 2+74 _ i_g+2745 . 8 
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